Convective gravity wave decoupling method based on physical information neural network and server

Abstract

The application relates to the technical field of computer data analysis, and provides a convection gravity wave decoupling method based on a physical information neural network and a server. The method performs joint decoupling analysis on a convection parameter field, a gravity wave activity field and an ozone concentration field of a target atmospheric region, separates the influence of dynamic transportation process and chemical process on the ozone distribution by using a neural network embedded with physical constraints, and determines the relative dominant role of different physical processes under different phases of the intra-seasonal oscillation by time sequence correlation analysis and causal inference, so as to provide an accurate basis for understanding the ozone exchange mechanism near the tropopause.

Description

Technical Field

[0001] This application belongs to the field of computer data analysis technology, specifically relating to a convective gravity wave decoupling method and server based on physical information neural networks. Background Technology

[0002] The distribution and variation of atmospheric ozone near the tropopause are jointly regulated by dynamic transport and chemical processes. Dynamic transport processes include horizontal advection transport caused by large-scale circulation, vertical upward transport induced by convective activity, and local mixing and exchange caused by gravity wave breakup. Chemical processes involve complex reaction pathways such as the transformation of ozone precursors generated by photochemical reactions, heterogeneous chemical depletion reactions, and the downward transport of stratospheric ozone to the troposphere. Under low-frequency variability conditions such as intraseasonal oscillations, the relative contributions of these two types of processes change significantly with phase evolution. Accurately quantifying the intensity of their influence on local ozone concentration variations is a key scientific issue in understanding the tropopause ozone budget balance.

[0003] Existing techniques typically rely on offline chemical transport models or statistical regression methods based on observational data to assess the dynamic and chemical contributions separately. The former is computationally expensive and highly sensitive to subgrid physical process parameterization schemes, while the latter struggles to distinguish between causal relationships and statistical correlations. Furthermore, existing methods lack effective means to simultaneously decouple the contributions of dynamic transport and chemical processes under unified physical constraints, leading to significant uncertainties in the attribution results. Summary of the Invention

[0004] This application provides a convective gravity wave decoupling method and server based on physical information neural networks.

[0005] This application provides a convective gravity wave decoupling method based on a physical information neural network, applied to a convective gravity wave decoupling server. The method includes: Acquire the set of convective parameter fields, gravity wave activity fields, and ozone concentration fields for the target atmospheric region; The set of convection parameters and the set of gravity wave activity fields are input into a preset physical information network for spatiotemporal evolution physical constraint processing, and the decoupled power transport contribution field set and the decoupled chemical process contribution field set are output. The spatiotemporal evolution physical constraint processing uses the Navier-Tex equation constraint term and the gravity wave linear theory constraint term as the physical regularization restriction condition in the network training process. The decoupled power transmission contribution field set and the ozone concentration field set are subjected to a first time series correlation analysis to extract a first time series response feature set between the power transmission contribution time series and the ozone concentration time series. The decoupled chemical process contribution field set and the ozone concentration field set are subjected to a second time series correlation analysis to extract a second time series response feature set between the chemical process contribution time series and the ozone concentration time series. The first time-series response feature set and the second time-series response feature set are input into a preset causal inference model to perform dynamic causal intensity inference processing to obtain a dynamic causal intensity distribution set, and the second time-series response feature set is then processed by the preset causal inference model to perform chemical causal intensity inference processing to obtain a chemical causal intensity distribution set; Relative dominance discrimination information is generated based on the set of dynamic causal intensity distributions and the set of chemical causal intensity distributions.

[0006] This application provides a convective gravity wave decoupling server, which includes a processor and a memory, wherein the memory stores a computer program, and when the computer program is executed by the processor, the processor performs the steps of the above method.

[0007] This application provides a computer-readable storage medium including a computer program. When the computer program is run on a convective gravity wave decoupling server, the computer program is used to cause the convective gravity wave decoupling server to perform the steps of the above-described method.

[0008] In this application, by embedding the constraints of the Navier-Tex equations and the linear theory of gravity waves into a pre-defined physical information network as physical regularization constraints, the spatiotemporal evolution physical constraints on the convective parameter field set and the gravity wave activity field set are processed. This achieves the decoupling and separation of the power transport process and the chemical process under the guidance of physical laws, ensuring that the output decoupled power transport contribution field set and decoupled chemical process contribution field set possess physical consistency and interpretability. Furthermore, by performing time-series correlation analysis on the decoupled contribution field set and the ozone concentration field set respectively, the time-series response feature sets between the power transport contribution time series and the chemical process contribution time series and the ozone concentration time series are extracted, effectively capturing the leading and lagging response relationships and attenuation characteristics of different physical processes to ozone concentration changes in the time dimension. Next, the two types of time-series response feature sets are input into a pre-defined causal inference model for causal intensity inference processing, obtaining the power causal intensity distribution set and the chemical causal intensity distribution set respectively, thereby quantitatively characterizing the spatial distribution of the causal influence intensity of power transport and chemical processes on ozone concentration changes. Finally, based on the two sets of causal intensity distribution sets, relative dominant role discrimination information was generated, which enabled an objective determination of the relative importance of the contribution of power transport and chemical process at different spatial locations, and improved the accuracy and reliability of ozone budget attribution analysis near the tropopause. Attached Figure Description

[0009] Figure 1 This is a flowchart illustrating a convective gravity wave decoupling method based on a physical information neural network, provided in an embodiment of this application.

[0010] Figure 2 This is a schematic diagram of the structure of a convective gravity wave decoupling server provided in an embodiment of this application.

[0011] Figure 3 This is a functional block diagram of a convective gravity wave decoupling server provided in an embodiment of this application. Detailed Implementation

[0012] To make the objectives, technical solutions, and advantages of the embodiments of this application clearer, the technical solutions of this application will be clearly and completely described below with reference to the accompanying drawings of the embodiments of this application. Obviously, the described embodiments are only some embodiments of the technical solutions of this application, and not all embodiments. Based on the embodiments recorded in this application, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the technical solutions of this application.

[0013] See Figure 1This is a convective gravity wave decoupling method based on physical information neural network provided in the embodiments of this application. This method can be applied to a convective gravity wave decoupling server. The specific process is as follows: steps 110-150.

[0014] This application provides a convective-gravity wave decoupling method based on a physical information neural network. This method performs joint decoupling analysis on the convective parameter field, gravity wave activity field, and ozone concentration field of the target atmospheric region. It uses a neural network with embedded physical constraints to separate the influence of dynamic transport processes and chemical processes on ozone distribution. Furthermore, it determines the relative dominant role of different physical processes in different phases of seasonal oscillations through time-series correlation analysis and causal inference, thereby providing an accurate basis for understanding the ozone exchange mechanism near the tropopause.

[0015] This application uses the decoupling analysis of convective gravity waves within the seasonal oscillation cycle of the atmospheric environment near the tropopause over a target atmospheric region as a detailed exemplary application scenario throughout the text. The target atmospheric region is selected as a three-dimensional spatial region covered by a latitude and longitude grid, vertically layered according to isobaric surfaces from the ground to a predetermined kilometer altitude, with a time span covering a complete seasonal oscillation life cycle. It should be understood that this example is only used to illustrate the technical implementation logic of this method and does not constitute a limitation on the scope of protection of this application.

[0016] Step 110: Obtain the set of convection parameters, the set of gravity wave activity, and the set of ozone concentration in the target atmospheric region.

[0017] The acquired set of convective parameter fields includes the convective effective potential energy values ​​of each spatial grid location in the target atmospheric region at different time steps. The convective effective potential energy is obtained by vertical integration of the work done by positive buoyancy during the adiabatic ascent of air parcels at each grid point.

[0018] The gravity wave activity field set includes the gravity wave flux divergence distribution at different time steps for each spatial grid location within the target atmospheric region. The gravity wave flux is calculated based on the cospectral distribution of zonal and vertical wind disturbances near the tropopause in wavenumber space, and a scalar field characterizing the upward propagation energy flux of gravity waves is obtained through averaging in the upper troposphere.

[0019] The ozone concentration field ensemble includes the ozone volume mixing ratio values ​​at different time steps for each spatial grid location within the target atmospheric region. It is obtained by reading the ozone mass mixing ratio data at the corresponding isobaric surface level in the reanalysis data and converting it into volume mixing ratio units.

[0020] The three field sets described above have the same number of time steps in the time dimension, denoted as T. They also have the same horizontal grid resolution in the spatial dimension, with the number of grids in the longitude direction denoted as Nx, the number of grids in the latitude direction as Ny, and the number of vertical layers denoted as Nz. All field sets are interpolated to a unified pressure layer coordinate system, with pressure layers distributed at varying intervals from 1000 hPa to 10 hPa. The acquired field sets are organized and stored as four-dimensional arrays, with each dimension corresponding to the time step index, vertical layer index, latitude index, and longitude index, respectively. The total data volume of a single field set is T × Nz × Ny × Nx floating-point values.

[0021] Step 120: Input the set of convection parameter fields and the set of gravity wave activity fields into the preset physical information network for spatiotemporal evolution physical constraint processing, and output the set of decoupled power transport contribution fields and the set of decoupled chemical process contribution fields. The spatiotemporal evolution physical constraint processing uses the Navier-Tex equation constraint terms and the gravity wave linear theory constraint terms as physical regularization constraints in the network training process.

[0022] Step 121: Extract convective effective potential energy distribution units at multiple spatial grid locations from the convective parameter field set, and extract gravity wave flux units at multiple spatial grid locations from the gravity wave activity field set.

[0023] The four-dimensional field ensemble is traversed along the spatial dimension at each time step. For each horizontal grid point, the convective effective potential energy values ​​of each pressure layer along its vertical profile are extracted, forming the convective effective potential energy distribution unit for that grid point at that time step. This distribution unit is a vector of length Nz. Simultaneously, for the same grid point, the gravity wave flux value at the same time step is extracted. This value has been weighted and averaged in the vertical direction to represent the upward propagation of gravity wave energy flux intensity at that grid point, also forming a vector of length Nz. The two vectors correspond one-to-one in the vertical layer dimension.

[0024] Step 122: Perform spatial grid coordinate alignment processing on the convective effective potential energy distribution unit and the gravitational fluctuation flux unit to obtain the aligned convective effective potential energy distribution and the aligned gravitational fluctuation flux distribution.

[0025] Since the convective parameter field and the gravity wave activity field may originate from different data assimilation systems or have different original grid resolutions, they need to be mapped to a preset standard grid coordinate system. The alignment process uses a bilinear interpolation algorithm, with the grid coordinates of the convective parameter field as the reference grid, and the gravity wave flux cells are interpolated to the grid point positions of the reference grid through latitude and longitude coordinate matching.

[0026] For each interpolation point on the reference grid, its latitude and longitude coordinates are obtained, and the four nearest neighbor source grid points surrounding these coordinates are located in the source grid coordinate system. Let the normalized relative position coordinates of the interpolation point in the source grid coordinate system be dx and dy, and the values ​​of the four nearest neighbor source grid points be V00, V10, V01, and V11, respectively. The calculation logic for the interpolation result V is as follows: First, perform two linear interpolations in the longitude direction to obtain the upper median value Vtop and the lower median value Vbot, where Vtop = V00 × (1 - dx) + V10 × dx and Vbot = V01 × (1 - dx) + V11 × dx. Then, perform one linear interpolation in the latitude direction to obtain V = Vtop × (1 - dy) + Vbot × dy. The above two-dimensional interpolation process is performed independently for each vertical layer of the gravity fluctuation flux unit. During interpolation, source grid data points located outside the reference grid are discarded, while locations covered by source grid data points within the reference grid are filled with nearest-neighbor valid values. After this processing, the two field sets have completely consistent spatial grid numbers and grid point coordinate correspondences at each time step.

[0027] Step 123: The aligned convective effective potential energy distribution and the aligned gravitational fluctuation flux distribution are sequentially combined along the time axis to generate a joint spatiotemporal grid data set.

[0028] The two aligned field sets are staggered and stacked along the time dimension. For each time step, the convective effective potential energy distribution is used as the first data layer, and the gravitational wave flux distribution is used as the second data layer. The two data layers correspond one-to-one in the spatial grid position, forming a dual-channel three-dimensional spatial data block. Each data block has a dimension of Nz×Ny×Nx×2. The dual-channel three-dimensional spatial data blocks of all time steps are connected end to end in chronological order to form a five-dimensional joint spatiotemporal grid data set with a dimension of T×Nz×Ny×Nx×2. This data set contains both thermodynamic parameters characterizing convection intensity and dynamic parameters characterizing gravitational wave activity, and the two are precisely aligned in the spatiotemporal coordinates.

[0029] Step 124: Call the shared feature extraction layer in the preset physical information network to perform spatial dimension downsampling and temporal dimension convolution processing on the joint spatiotemporal grid data set to obtain a spatiotemporal coupled potential feature representation set. Call the power transmission mapping branch in the preset physical information network to perform grid-by-grid feature transformation processing on the spatiotemporal coupled potential feature representation set to generate a preliminary distribution set of power transmission contributions. Call the chemical process mapping branch in the preset physical information network to perform grid-by-grid feature transformation processing on the spatiotemporal coupled potential feature representation set to generate a preliminary distribution set of chemical process contributions.

[0030] The pre-defined physical information network adopts an encoder-dual-branch decoder architecture. The shared feature extraction layer consists of three cascaded 3D convolutional blocks.

[0031] The first convolutional block contains one 3D convolutional layer, one batch normalization layer, and one linear unit activation layer with leakage correction. The kernel size of the 3D convolutional layer is set to 3×3×3, the convolution stride is set to 2 in the spatial dimension and 1 in the temporal dimension, and the number of output feature channels is set to 32.

[0032] The structure of the second convolutional block is the same as that of the first convolutional block. The kernel size of the three-dimensional convolutional layer is set to 3×3×3, the convolution stride is set to 2 in the spatial dimension and 1 in the temporal dimension, the number of input feature channels is 32, and the number of output feature channels is 64.

[0033] The structure of the third convolutional block is the same as the previous convolutional blocks. The kernel size of the three-dimensional convolutional layer is set to 3×3×3, the convolution stride is set to 1 in the spatial dimension and 1 in the temporal dimension, the number of input feature channels is 64, and the number of output feature channels is set to 128.

[0034] The data flow process of the joint spatiotemporal grid dataset through the above three convolutional blocks is as follows: The input five-dimensional data first enters the three-dimensional convolutional layer of the first convolutional block. The convolutional kernel slides with a stride of 2 in the spatial dimension of the input data and with a stride of 1 in the temporal dimension. At each slide, the sum of the element-wise products of the convolutional kernel weight and the corresponding input data segment is calculated, and a scalar value is output after adding the bias term. The output values ​​of all sliding positions are combined to form an intermediate feature map. The intermediate feature map is processed by a batch normalization layer to normalize the channel direction mean and variance. The normalization result is processed by a linear unit activation layer with leakage correction for nonlinear transformation. The negative value part of the output is scaled according to a preset slope coefficient, and the positive value part is passed directly.

[0035] The first convolutional block outputs a feature representation with a dimension of T×floor(Nz / 2)×floor(Ny / 2)×floor(Nx / 2)×32. This feature representation enters the second convolutional block and undergoes the same convolution, normalization, and activation operations, further reducing the spatial resolution to one-quarter of the original resolution. The output feature representation has a dimension of T×floor(Nz / 4)×floor(Ny / 4)×floor(Nx / 4)×64. This feature representation enters the third convolutional block. Since the spatial stride is set to 1, the spatial resolution remains unchanged, and the temporal dimension length remains unchanged. The output is a set of spatiotemporally coupled latent feature representations with a dimension of T×floor(Nz / 4)×floor(Ny / 4)×floor(Nx / 4)×128.

[0036] The power transmission mapping branch consists of two pointwise convolutional layers and one trilinear upsampling layer. The kernel size of the pointwise convolutional layers is 1×1×1, used to perform independent linear transformations on the feature vectors at each spatial grid location without introducing spatial neighborhood information mixing. The spatiotemporally coupled latent feature representation set first enters the first pointwise convolutional layer, which has 128 input channels and 64 output channels. Each output channel corresponds to a 1×1×1 convolutional kernel, which performs a weighted summation of the input feature vectors along the channel direction, outputting an intermediate representation with a dimension of T×floor(Nz / 4)×floor(Ny / 4)×floor(Nx / 4)×64. This intermediate representation enters the second pointwise convolutional layer, which has 64 input channels and 1 output channel, outputting an initial contribution distribution with a dimension of T×floor(Nz / 4)×floor(Ny / 4)×floor(Nx / 4)×1.

[0037] The trilinear upsampling layer interpolates the initial contribution distribution across three spatial dimensions to restore it to the same grid size as the original input field set. The interpolation method used is trilinear interpolation, performing linear interpolation calculations sequentially in the vertical, latitude, and longitude directions. The output after upsampling is a preliminary distribution set of the power transport contribution, with dimensions T×Nz×Ny×Nx×1, and its value represents the local ozone variation tendency term caused by the power transport process.

[0038] The chemical process mapping branch has the exact same structural configuration as the dynamic transport mapping branch, consisting of two pointwise convolutional layers and one trilinear upsampling layer. The first pointwise convolutional layer has 128 input channels and 64 output channels, and the second pointwise convolutional layer has 64 input channels and 1 output channel. The upsampling layer uses trilinear interpolation. The convolutional kernel weight parameters within the chemical process mapping branch are independent of those in the dynamic transport mapping branch and are learned separately during network training. The preliminary distribution set of chemical process contributions output by the chemical process mapping branch also has dimensions of T×Nz×Ny×Nx×1, and its value represents the ozone localization tendency term caused by chemical processes such as photochemical reactions.

[0039] Step 125: Construct a physical rule compliance evaluation function based on the preliminary distribution set of power transport contribution and the preliminary distribution set of chemical process contribution. The physical rule compliance evaluation function includes a compliance evaluation sub-function of the Navier-Tex equation and a compliance evaluation sub-function of the linear theory of gravity waves.

[0040] Step 1251: Separate the velocity field divergence component distribution unit and the momentum field advection component distribution unit from the initial distribution set of power transport contribution.

[0041] The output values ​​of the preliminary distribution set of the power transport contribution include information on the transport effect of the horizontal wind field on the ozone field. Based on the dynamic expansion of the ozone continuity equation, the power transport contribution can be decomposed into two parts: a velocity field divergence term and an advection transport term. The velocity field divergence component is proportional to the convergence and divergence of the horizontal wind field, while the momentum field advection component is proportional to the migration of the ozone background gradient under the influence of the wind field. Separation processing is achieved by linearly decoupling the single-channel results output from the power transport mapping branch from the pre-calculated horizontal wind field divergence field and ozone background gradient field grid by grid.

[0042] For each spatial grid location, an independent variable matrix is ​​constructed, consisting of the horizontal wind divergence values ​​and ozone background gradient values ​​of that grid and its adjacent grids. The preliminary distribution of the contribution of power transport is used as the dependent variable, and the decoupling coefficients are independently solved at each time step using a least-squares fitting method. The decoupling coefficients include two scalar values: divergence coefficient and advection coefficient. The divergence coefficient multiplied by the horizontal wind divergence field yields the velocity field divergence component distribution unit, and the advection coefficient multiplied by the projection of the ozone background gradient field onto the horizontal wind direction yields the momentum field advection component distribution unit.

[0043] Step 1252: Analyze the divergence deviation distribution between the velocity field divergence component distribution unit and the incompressible fluid continuity condition to obtain the continuity constraint compliance component of the Navier-Tox equation. Analyze the momentum balance deviation distribution between the momentum field advection component distribution unit and the gravity fluctuation source-sink term to obtain the momentum constraint compliance component of the Navier-Tox equation.

[0044] For the continuity constraint compliance component, at each time step, the absolute value of the velocity field divergence component at each spatial grid location is taken. The continuity condition for incompressible fluids requires that the velocity field divergence is always equal to 0; therefore, the smaller the absolute value of the divergence component, the closer it is to the zero divergence condition of incompressible fluids. The absolute values ​​of divergence at all spatial grid locations are averaged in an area-weighted manner on a horizontal plane, with the latitudinal weights using the cosine of the latitude corresponding to the grid. The resulting weighted average is the continuity constraint compliance component value for that time step.

[0045] For the momentum constraint compliance component, at each time step, the values ​​of the momentum field advection component and the gravity wave flux divergence field at the corresponding grid location are calculated point-by-point. The gravity wave flux divergence field is obtained by taking the central difference of the gravity wave flux cells in the vertical direction within the gravity wave activity field set. For each spatial grid location, the difference between the momentum field advection component value and the gravity wave flux divergence value is calculated. The sum of the squares of the differences for all grid locations is divided by the total number of effective grids, and the square root is taken to obtain the root mean square value of the momentum balance deviation. The smaller the root mean square value, the higher the degree of matching between the momentum field advection component and the gravity wave source and sink terms. The reciprocal of this root mean square value is normalized and used as the momentum constraint compliance component value.

[0046] Step 1253: Separate the vertical wavenumber characteristic distribution and the horizontal wavenumber characteristic distribution of gravity waves from the gravity wave flux unit of the gravity wave activity field set, and extract the background flow field adjustment component distribution caused by the power transport process from the preliminary distribution set of power transport contribution.

[0047] The vertical wavenumber characteristic distribution of gravity waves is obtained by performing Fourier transform decomposition on the vertical profile of the gravity wave flux unit. For each horizontal grid location and each time step, the vertical profile vector at that location is extracted, and a Fast Fourier Transform algorithm is applied to this vector to obtain the complex amplitude spectrum in the vertical wavenumber domain. The wavenumber component with the largest amplitude modulus is searched in the wavenumber spectrum; the wavenumber value corresponding to this wavenumber is the characteristic value of the vertical wavenumber of gravity waves at that grid location. The above operation is performed for all horizontal grid locations and all time steps to obtain the characteristic distribution of the vertical wavenumber of gravity waves with dimensions T×Ny×Nx.

[0048] The horizontal wavenumber characteristic distribution of gravity waves is obtained by performing a two-dimensional Fourier transform on the horizontal distribution field of gravity wave flux at the same time step. For each time step and each vertical layer, a two-dimensional scalar field on the horizontal plane of that layer is extracted, and a two-dimensional fast Fourier transform algorithm is applied to obtain the two-dimensional complex amplitude spectra in the zonal and meridional wavenumber domains. The vector composite modulus of the zonal and meridional wavenumber components is calculated, i.e., the square root of the sum of the squares of the zonal and meridional wavenumbers, and is taken as the horizontal wavenumber modulus of that wavenumber point. The wavenumber point with the largest amplitude modulus is searched in the two-dimensional wavenumber spectrum; the horizontal wavenumber modulus corresponding to this position is the characteristic value of the horizontal wavenumber of gravity waves at that time step and vertical layer. The above operation is performed for all time steps and all vertical layers to obtain the characteristic distribution of the horizontal wavenumber of gravity waves.

[0049] The background flow field adjustment component distribution characterizes the feedback correction effect of the power transport process on the background horizontal wind field. The preliminary distribution set of power transport contributions is input into a pre-trained linear regression model to obtain the background flow field adjustment component distribution. This linear regression model is trained using the ozone tendency time series over a historical period as input features and the wind field tendency time series as the prediction target. The model weight matrix has a dimension of 1×M, where M is the number of historical time steps used for prediction. During training, the least squares method is used to solve the weight matrix to minimize the mean square error between the model's predicted wind field tendency and the actual observed wind field tendency. After training, the time series of the preliminary distribution set of power transport contributions is input into the linear regression model, and the wind field tendency estimate output by the model is the background flow field adjustment component distribution.

[0050] Step 1254: Based on the vertical wavenumber characteristic distribution of gravity waves, the horizontal wavenumber characteristic distribution of gravity waves, and the background flow field adjustment component distribution, analyze the dispersion relation response deviation distribution and the refraction deviation distribution of energy propagation direction of gravity waves to obtain the dispersion constraint compliance component and the refraction constraint compliance component of the linear theory of gravity waves.

[0051] In the linear theory of gravity waves, the frequency and wavenumber of gravity waves satisfy a dispersion relation, i.e., the square of the natural frequency equals the square of the buoyancy frequency multiplied by the ratio of the square of the horizontal wavenumber to the square of the total wavenumber. The square of the total wavenumber is the sum of the squares of the horizontal and vertical wavenumbers, and the square of the vertical wavenumber equals the vertical wavenumber multiplied by 2π divided by the vertical wavelength, squared. For each spatial grid location and each time step, the theoretical natural frequency is calculated by substituting the extracted vertical and horizontal wavenumber characteristic values, combined with the square of the buoyancy frequency at that grid location, into the dispersion relation expression. The absolute value of the difference between this theoretical natural frequency and the actual observed gravity wave frequency is the dispersion relation response deviation. The actual observed gravity wave frequency is obtained by performing wavelet transform analysis on the time series of the gravity wave activity field set. The average value of the dispersion relation response deviation for all grid locations and time steps is calculated to obtain the dispersion constraint compliance component value.

[0052] The refraction constraint compliance component is calculated based on the degree of consistency between the gravity wave group velocity direction and the background flow field adjustment direction. The group velocity direction is determined by the partial derivative of the natural frequency with respect to the wavenumber vector; specifically, the horizontal component of the group velocity is equal to the horizontal group velocity of the gravity wave, and its direction is parallel to the horizontal wavenumber vector direction. The background flow field adjustment direction is determined by the gradient direction of the background flow field adjustment component distribution, which is obtained by calculating the first-order central difference of the background flow field adjustment component in the zonal and meridional directions. For each spatial grid location, the dot product of the group velocity unit vector and the background flow field adjustment direction unit vector is calculated. This dot product value is equal to the cosine of the angle between the two direction vectors. The cosine value is 1 when the two directions are completely consistent and -1 when they are completely opposite. The cosine value is averaged across all spatial grid locations. The closer the average value is to 1, the higher the consistency between the gravity wave energy propagation direction and the background flow field adjustment direction. The normalized result of this average value constitutes the refraction constraint compliance component value.

[0053] Step 1255: Spatial averaging and fusion processing of the continuity constraint compliance component and the momentum constraint compliance component is performed to generate the compliance evaluation sub-function of the Navier-Tox equation.

[0054] Spatial averaging fusion processing employs a weighted averaging method, with weight coefficients adjusted for Coriolis parameters based on the latitudinal location of each grid point. The Coriolis parameter f = 2 × Ω × sinφ, where Ω is the Earth's rotation angular rate and φ is the latitudinal radian value. The weight of the continuity constraint compliance component is set to 1, and the weight of the momentum constraint compliance component is set to the square of f divided by the square of the maximum value of f, where the maximum value of f is the maximum absolute value of the Coriolis parameter within the latitudinal range of the target atmospheric region. For each time step, the continuity constraint compliance component at each spatial grid location is multiplied by its weight of 1, and the momentum constraint compliance component is multiplied by its corresponding latitudinal weight. The sum of these values ​​is divided by the total sum of weights to obtain the fusion scalar value for that time step. This scalar value is the output value of the Navier-Tox equation compliance evaluation subfunction for that time step. The smaller this value, the more reasonable the decoupling result output by the preset physical information network is at the Navier-Tox equation constraint level.

[0055] Step 1256: Spatial averaging and fusion of the dispersion constraint compliance component and the refraction constraint compliance component are performed to generate a compliance evaluation sub-function for the linear theory of gravity waves.

[0056] Spatial averaging fusion employs a weighted averaging method. The weights of the dispersion constraint compliance components are adjusted based on the vertical wavenumber magnitude. For each spatial grid location, the ratio of the characteristic vertical wavenumber value at that location to a preset reference vertical wavenumber is calculated. This ratio is used as the weight coefficient for the dispersion constraint compliance component, with higher weights assigned to regions with larger vertical wavenumbers to reflect the greater sensitivity of short-wave gravity waves to dispersion relations. The weight of the refractive constraint compliance component is set to 1. For each time step, the dispersion constraint compliance component at each spatial grid location is multiplied by its corresponding vertical wavenumber weight, and the refractive constraint compliance component is multiplied by its weight of 1. The sum of these values ​​is divided by the total weight sum to obtain the fusion scalar value for that time step. This scalar value is the output value of the gravity wave linear theory compliance evaluation sub-function for that time step.

[0057] Step 1257: Perform a weighted summation of the conformity evaluation sub-function of the Navier-Tox equation and the conformity evaluation sub-function of the linear theory of gravity waves to generate the output value of the physical rule conformity evaluation function.

[0058] The weighting coefficients for both sub-functions are set to 0.5, indicating that fluid dynamics constraints and wave dynamics constraints have equally important regularization status during network training. For each time step, the output value of the physical rule compliance evaluation function is equal to 0.5 multiplied by the output value of the Navier-Tox equation compliance evaluation sub-function plus 0.5 multiplied by the output value of the gravity wave linear theory compliance evaluation sub-function. The average of the output values ​​over all time steps yields the physical rule compliance evaluation function output value for the entire time series. The smaller this value, the more reasonable the decoupling result of the preset physical information network output is at the physical rule level.

[0059] Step 126: Construct a joint optimization objective function based on the output value of the physical rule compliance evaluation function and the data reconstruction difference value, and iteratively adjust the internal connection weights of the preset physical information network according to the joint optimization objective function.

[0060] The data reconstruction difference value (MSE) is the sum of the preliminary distribution sets of power transport contributions and chemical process contributions output by the preset physical information network, and the mean square error between this sum and the observed ozone concentration temporal trend. The observed ozone concentration temporal trend is obtained by performing a central difference operation on the ozone concentration field set along the time dimension. For each spatial grid location and each time step, the ozone concentration temporal trend is equal to the difference between the ozone concentration value of the next time step and the ozone concentration value of the previous time step, divided by twice the time step size. Let the element value of the preliminary distribution set of power transport contributions be Vadv, the element value of the preliminary distribution set of chemical process contributions be Vchem, and the observed value of the ozone concentration temporal trend be Tobs. Then, the error for a single grid location is Vadv + Vchem - Tobs. The sum of the squared errors for all spatial grid locations and all time steps, divided by the total number of valid data points, yields the data reconstruction difference value (MSE).

[0061] The joint optimization objective function is constructed as Ltotal = MSE + λ × Lphy, where Lphy is the output value of the physical rule compliance evaluation function, and λ is the balance coefficient, empirically set to 0.1. Iterative adjustment employs an adaptive moment estimation optimization algorithm, which maintains the first and second moment estimates for each network weight parameter and calculates the weight update amount using the bias-corrected moment estimates. The initial learning rate is set to 0.001, the first moment decay rate to 0.9, the second moment decay rate to 0.999, and the numerical stability constant to 1 × 10⁻⁶. -8 In each iteration, a batch of the joint spatiotemporal grid dataset is input into a predefined physical information network. Forward propagation calculates the outputs of each network layer, and the loss value is calculated based on the joint optimization objective function. Backpropagation calculates the gradient of the loss value with respect to the weights of each network layer, and the adaptive moment estimation optimizer updates the weight values ​​based on the gradient direction. Every 100 iterations, the learning rate decays to 0.9 times the previous stage. Training continues until the relative decrease in the joint optimization objective function is less than 0.001 within 20 consecutive iterations, at which point the iterative adjustment process is considered to have reached a preset stable state.

[0062] Step 127: When the iteration adjustment process reaches a preset stable state, the preliminary distribution set of power transmission contribution output by the power transmission mapping branch is output as the decoupled power transmission contribution field set, and the preliminary distribution set of chemical process contribution output by the chemical process mapping branch is output as the decoupled chemical process contribution field set.

[0063] Once a preset stable state is reached, the network weight parameters are no longer updated. The joint spatiotemporal grid data set of the entire time series is input step-by-step into the trained preset physical information network. The network performs forward propagation computation, and the shared feature extraction layer performs spatial downsampling and feature extraction on the input data. The dynamic transport mapping branch and the chemical process mapping branch transform and upsample the extracted latent features, respectively. The outputs of the two branches are collected at all time steps and concatenated in the original time order to obtain the complete decoupled dynamic transport contribution field set and the decoupled chemical process contribution field set. The spatial extent, grid resolution, and number of time steps of the two sets are consistent with the input convection parameter field set and gravity wave activity field set, and the dimensions are T×Nz×Ny×Nx×1.

[0064] Step 130: Perform a first time-series correlation analysis on the decoupled power transmission contribution field set and the ozone concentration field set to extract a first time-series response feature set between the power transmission contribution time series and the ozone concentration time series; and perform a second time-series correlation analysis on the decoupled chemical process contribution field set and the ozone concentration field set to extract a second time-series response feature set between the chemical process contribution time series and the ozone concentration time series.

[0065] Step 131: Extract time series units of power transport contribution intensity arranged along the time direction from the decoupled power transport contribution field set for each spatial grid, and extract time series units of ozone concentration value arranged along the time direction from the ozone concentration field set for the same spatial grid. Perform synchronous sliding time window division processing on the time series units of power transport contribution intensity and ozone concentration value, to obtain multiple sets of time window segments corresponding to power transport contribution window segments and ozone concentration window segments.

[0066] For each spatial grid location within the target atmospheric region, in a T×Nz×Ny×Nx×1 dimensional array of the decoupled dynamic transport contribution field set, with fixed vertical layer indices, latitude indices, and longitude indices, a one-dimensional vector of length T is extracted along the time dimension. This vector represents the time series unit of the dynamic transport contribution intensity at that spatial grid location. Using the same index location, a one-dimensional vector of length T is extracted from the ozone concentration field set to obtain the ozone concentration value time series unit.

[0067] The sliding time window partitioning process uses a window length of W time steps and a sliding step size of S time steps. For a time series unit of length T, starting from time step index 1, the subsequence from the 1st time step to the Wth time step is extracted using a window length of W as the first group of power transport contribution window segments and the first group of ozone concentration window segments. Then, the window slides forward S time steps, extracting the subsequence from the (1+S)th time step to the (W+S)th time step as the second group of window segments.

[0068] This process continues until the end of the window exceeds the total length of the time series. The number of window groups N obtained from the partitioning process is equal to floor((TW)÷S)+1. For time segments containing missing values ​​within the window, linear interpolation is used for imputation. The imputation formula is: the missing value is equal to the difference between the previous valid value and the next valid value minus the previous valid value, multiplied by the distance between the missing position and the previous valid position, divided by the total distance between the two valid positions.

[0069] Step 132: Perform hysteresis cross-correlation analysis on the power transport contribution window segment and ozone concentration window segment corresponding to each time window to generate the first hysteresis correlation coefficient variation curve at each spatial grid location. Identify the correlation coefficient peak points in the positive hysteresis region and the negative hysteresis region from the first hysteresis correlation coefficient variation curve, and extract the first lead response time parameter corresponding to the positive hysteresis peak point and the first hysteresis response time parameter corresponding to the negative hysteresis peak point.

[0070] The implementation of lag cross-correlation analysis is as follows: For a given number of lag steps L, where L is an integer ranging from -(W+1) to (W-1), the power transport contribution window sequence is shifted L time steps along the time axis. When L is greater than 0, the power transport contribution window is shifted backward, indicating that the change in power transport leads the change in ozone concentration by L time steps. When L is less than 0, the power transport contribution window is shifted forward, indicating that the change in power transport lags the change in ozone concentration by |L| time steps. The effective length of the two sequences in the overlapping time steps after shifting is W-|L|. The Pearson correlation coefficient between the two sequences in the overlapping portion after shifting is calculated. The correlation coefficient r is calculated as the product of the covariance of the two sequences in the overlapping portion and the standard deviation of the two sequences.

[0071] For each value of L, a correlation coefficient is calculated. Plotting L on the horizontal axis and the correlation coefficient on the vertical axis, connect the correlation coefficients corresponding to all L values ​​to form a curve. This curve represents the first lag correlation coefficient variation curve at that spatial grid location. Within the positive lag region where L is greater than 0, search for the location of the maximum correlation coefficient. The lag step number Lmax1 at this location, multiplied by the duration of a single time step Δt, is the first lead response time parameter, in time units. Within the negative lag region where L is less than 0, search for the location of the maximum correlation coefficient. The absolute value of the lag step number |Lmax2| at this location, multiplied by Δt, is the first lag response time parameter.

[0072] Step 133: Analyze the correlation coefficient decay trend on both sides of the peak point of the first lag correlation coefficient change curve as the absolute value of the lag increases, and extract the first correlation coefficient decay characteristic parameter that characterizes the rate of decay.

[0073] Centered on the peak point Lmax1 of the positive lag region, intervals of R lag steps are taken to the left and right, where R is the preset fitting radius. The correlation coefficient value sequence within this interval is fitted using an exponential function, with the fitting function form r(L) = rmax × exp(-K × |L - Lmax1|), where rmax is the peak correlation coefficient, |L - Lmax1| is the absolute value of the lag step offset from the peak point, and K is the decay coefficient to be fitted. During the fitting process, the natural logarithm of both sides of the fitting function is taken, transforming the exponential fitting into a linear fitting problem, i.e., ln(r(L)) = ln(rmax) - K × |L - Lmax1|. The value of K is solved using the least squares method; this K value is the first correlation coefficient decay characteristic parameter. A larger K value indicates that the correlation coefficient decays faster with increasing lag offset, the correlation's temporal persistence is weaker, and the impact of the power transmission process on ozone concentration is more concentrated in time.

[0074] Step 134: Combine the first lead response time parameter, the first lag response time parameter, and the first correlation coefficient decay characteristic parameter into the first time series response feature set corresponding to each spatial grid.

[0075] For each spatial grid location, the above three parameters are concatenated into a 3D feature vector in the order of the first lead response time parameter, the first lag response time parameter, and the first correlation coefficient decay characteristic parameter. The feature vectors of all Nz×Ny×Nx spatial grid locations within the target atmospheric region together constitute the first temporal response feature set. This set is a four-dimensional array in the spatial dimension, with each dimension being the vertical layer index, latitude index, longitude index, and feature index, respectively. The array size is Nz×Ny×Nx×3.

[0076] Step 135: Extract time series units of chemical process contribution intensity arranged along the time direction for each spatial grid from the decoupled chemical process contribution field set. Perform synchronous sliding time window division on the time series units of chemical process contribution intensity and ozone concentration value, to obtain multiple sets of time window segments corresponding to chemical process contribution window segments and ozone concentration window segments.

[0077] The extraction method for the time series units of chemical process contribution intensity is completely consistent with that for the extraction of the time series units of power transport contribution intensity in step 131. A one-dimensional vector of length T is extracted from the decoupled chemical process contribution field set using the same spatial index. The window length W and sliding step size S parameters for the sliding time window partitioning process are set exactly the same as in step 131 to ensure that the two sets of time series correlation analysis processes have comparable time analysis scales. After partitioning, N sets of chemical process contribution window segments and ozone concentration window segments are also obtained at each spatial grid location.

[0078] Step 136: Perform hysteresis cross-correlation analysis on the chemical process contribution window segment and ozone concentration window segment corresponding to each time window to generate the second hysteresis correlation coefficient variation curve at each spatial grid location. Identify the correlation coefficient peak points in the positive hysteresis region and the negative hysteresis region from the second hysteresis correlation coefficient variation curve, and extract the second lead response time parameter corresponding to the positive hysteresis peak point and the second hysteresis response time parameter corresponding to the negative hysteresis peak point.

[0079] The calculation process for lag cross-correlation analysis is completely consistent with step 132. For each lag step L, the Pearson correlation coefficient between the chemical process contribution window and the ozone concentration window is calculated, generating a second lag correlation coefficient variation curve. Within the positive lag region where L is greater than 0, the location of the maximum correlation coefficient point is searched; the lag step number corresponding to this location multiplied by Δt is the second lead response time parameter. Within the negative lag region where L is less than 0, the location of the maximum correlation coefficient point is searched; the absolute value of the lag step number corresponding to this location multiplied by Δt is the second lag response time parameter. The second lead response time parameter and the second lag response time parameter respectively characterize the typical time scales at which changes in the chemical process lead and lag behind changes in ozone concentration.

[0080] Step 137: Analyze the correlation coefficient decay trend on both sides of the peak point of the second lag correlation coefficient change curve as the absolute value of the lag increases, and extract the second correlation coefficient decay characteristic parameter that characterizes the rate of decay.

[0081] The method for extracting the attenuation characteristic parameter is consistent with the exponential function fitting method described in step 133. Taking the peak point of the positive lag region of the second lag correlation coefficient change curve as the center, the correlation coefficient value sequence within the fitting radius R is subjected to exponential function fitting. The attenuation coefficient K obtained by fitting is the attenuation characteristic parameter of the second correlation coefficient, which characterizes the time duration of the effect of the chemical process on ozone concentration.

[0082] Step 138: Combine the second lead response time parameter, the second lag response time parameter, and the second correlation coefficient decay characteristic parameter into a second time series response feature set corresponding to each spatial grid.

[0083] The combination method is the same as in step 134, where the second lead response time parameter, the second lag response time parameter, and the second correlation coefficient decay characteristic parameter are concatenated into a 3D feature vector. The feature vectors of all Nz×Ny×Nx spatial grids constitute the second time-series response feature set, with an array size of Nz×Ny×Nx×3.

[0084] Step 140: Input the first time-series response feature set and the second time-series response feature set into the preset causal inference model to perform dynamic causal intensity inference processing to obtain the dynamic causal intensity distribution set, and perform chemical causal intensity inference processing on the second time-series response feature set through the preset causal inference model to obtain the chemical causal intensity distribution set.

[0085] Step 141: The first response feature vector corresponding to each spatial grid in the first time-series response feature set is concatenated with the second response feature vector corresponding to the same spatial grid in the second time-series response feature set to generate a fused causal feature representation set of multiple spatial grids.

[0086] The first response feature vector is a 3-dimensional vector containing the first lead response time parameter, the first lag response time parameter, and the first correlation coefficient decay characteristic parameter. The second response feature vector is also a 3-dimensional vector containing the second lead response time parameter, the second lag response time parameter, and the second correlation coefficient decay characteristic parameter. The two sets of vectors are concatenated end-to-end along their feature dimensions, resulting in a 6-dimensional fused causal feature vector for each spatial grid position. All Nz×Ny×Nx fused causal feature vectors from these spatial grids are arranged in their original spatial positions, forming a fused causal feature representation set with an array size of Nz×Ny×Nx×6.

[0087] Step 142: Call the first causal test branch in the preset causal inference model to perform a causal significance assessment of the fused causal feature representation set, analyze the degree of improvement in the prediction ability of ozone concentration time series before and after the introduction of the power transmission contribution time series, and obtain the distribution of power transmission causal significance index.

[0088] Step 1421: Decompose the historical observation sequence of power transport contribution and the historical observation sequence of ozone concentration from the fused causal feature representation set. Based on the historical observation sequence of ozone concentration, construct a restricted prediction model that only relies on the historical information of ozone concentration itself. Perform forward extrapolation processing on the current ozone concentration through the restricted prediction model to obtain the restricted ozone concentration extrapolation sequence.

[0089] For each spatial grid location, the historical observation sequence of power transport contribution, Adv(t), is obtained from the decoupled power transport contribution field set, where t ranges from 1 to T. The historical observation sequence of ozone concentration, O3(t), is obtained from the ozone concentration field set. The constrained prediction model uses a vector autoregression model, and the model order p is automatically determined according to the Bayesian information criterion. The Bayesian information criterion is calculated as BIC = T × ln(σ). 2 ) + p × ln(T), where σ 2This represents the model residual variance. The model order is determined by iterating through p from 1 to the preset maximum order Pmax, calculating the BIC value for each p value, and selecting the p that minimizes the BIC value.

[0090] The constrained prediction model takes the form O3(t) = c + Σ(i=1top)φi×O3(ti) + ε(t), where c is a constant term, φi is the autoregressive coefficient, and ε(t) is the white noise residual term. The model coefficients c and φi are obtained by least squares fitting, with the fitting objective being to minimize the sum of squared residuals Σ(t=p+1toT)[O3(t) - c - Σ(i=1top)φi×O3(ti)] 2 After fitting, the historical ozone concentration observation sequence is slidably input into the trained restricted prediction model step by step. For each time step t, p historical values ​​from O3(tp) to O3(t-1) are input, and the model outputs the ozone concentration projection value for that time step. The projection values ​​of all time steps constitute the complete restricted ozone concentration projection sequence O3_restricted(t).

[0091] Step 1422: Obtain an index of the degree of difference between the restricted ozone concentration projection sequence and the actual ozone concentration observation sequence.

[0092] The root mean square error (RMSE) is used to calculate the degree of difference in the restricted ozone concentration projection. For a sequence length T, the RMSE-restricted index is calculated as follows: First, the square of the point-by-point difference between the restricted ozone concentration projection sequence O3_restricted(t) and the actual ozone concentration observation sequence O3(t) is calculated, i.e., [O3_restricted(t) - O3(t)]. 2 The summation of the squared differences over all valid time steps t is divided by the total number of valid time steps Tvalid, and the square root of the quotient is then taken. Tvalid is the number of time steps that the constrained prediction model can extrapolate, and is equal to Tp.

[0093] Step 1423: Based on the historical observation sequence of power transport contribution and the historical observation sequence of ozone concentration, construct an unconstrained prediction model that includes external inputs. Use the unconstrained prediction model to perform forward extrapolation on the current ozone concentration to obtain the unconstrained ozone concentration extrapolation sequence.

[0094] The unrestricted prediction model employs a vector autoregressive model with exogenous variables. The model inputs include the past p time steps of historical ozone concentration observations and the past q time steps of historical power transport contribution observations. The model order p and q are both determined using the Bayesian information criterion, which is calculated as BIC = T × ln(σ). 2) + (p+q)×ln(T). Iterate through p and q from 1 to Pmax and calculate the BIC value for each combination. Select the combination that minimizes the BIC value as the final model order.

[0095] The unconstrained prediction model takes the form O3(t) = c + Σ(i=1 to top)φi×O3(ti) + Σ(j=1 to q)ψj×Adv(tj) + ε(t), where c is a constant term, φi is the autoregressive coefficient, ψj is the exogenous variable coefficient, and ε(t) is the white noise residual term. The model coefficients are obtained by least squares fitting, with the fitting objective being to minimize the sum of squared residuals Σ(t=max(p,q)+1 to T)[O3(t)-c-Σ(i=1 to top)φi×O3(ti)-Σ(j=1 to q)ψj×Adv(tj)] 2 After fitting, the historical values ​​of ozone concentration and power transport contribution are input step by step to obtain the unrestricted ozone concentration projection sequence O3_unrestricted(t).

[0096] Step 1424: Obtain an index of the degree of difference between the unrestricted ozone concentration projection sequence and the actual ozone concentration observation sequence.

[0097] The calculation method is the same as in step 1422. The calculation method for the unrestricted inference difference index RMSEunrestricted is as follows: calculate the sum of squares of the point-by-point differences between O3_unrestricted(t) and O3(t), divide by the total number of effective time steps, and then take the square root. The total number of effective time steps is T-max(p, q).

[0098] Step 1425: Construct a causal significance evaluation function for power transport based on the relative change magnitude between the constrained and unconstrained inference degree indices, and obtain the causal significance evaluation value for power transport for each spatial grid.

[0099] The causal significance evaluation function for power transmission is defined as follows: Fscore = (RMSErestricted - RMSEunrestricted) ÷ RMSErestricted. This evaluation value represents the relative reduction in prediction error after incorporating power transport contribution information. An Fscore greater than 0 indicates improved prediction accuracy after incorporating power transport contribution information; a larger Fscore indicates a stronger explanatory power of power transport contribution to ozone concentration changes. If RMSEunrestricted is greater than or equal to RMSErestricted, then Fscore is assigned a value of 0.

[0100] Step 1426: Compare the power transmission causal significance assessment value with the preset significance discrimination conditions, and determine the distribution of power transmission causal significance index based on the comparison results.

[0101] The preset significance criteria are as follows: grid positions with Fscore greater than Thigh are marked as having significant causal relationship in power transmission, and the causal significance index for power transmission is assigned a value of 1. Grid positions with Fscore less than or equal to Thigh but greater than Tlow are marked as having weakly significant causal relationship in power transmission, and the index is assigned a value of 0.5. Grid positions with Fscore less than or equal to Tlow are marked as having insignificant causal relationship in power transmission, and the index is assigned a value of 0. Thigh and Tlow are preset threshold parameters. The index assignment results for all Nz×Ny×Nx spatial grids constitute the causal significance index distribution for power transmission, with an array size of Nz×Ny×Nx.

[0102] Step 143: Call the second causal test branch in the preset causal inference model to perform chemical process causal significance assessment on the fused causal feature representation set, analyze the degree of improvement in the prediction ability of ozone concentration time series before and after the introduction of chemical process contribution time series, and obtain the distribution of chemical process causal significance index.

[0103] The specific implementation methods for steps 1431 to 1437 are the same as those for steps 1421 to 1426, except that the input exogenous variable is replaced by the historical observation sequence Adv(t) of the contribution of dynamic transport with the historical observation sequence Chem(t) of the contribution of chemical processes. The expression for the causal significance assessment function of the chemical process is the same:

[0104] Gscore=(RMSErestricted-RMSEunrestricted)÷RMSErestricted.

[0105] The threshold parameters Thigh and Tlow for the preset significance discrimination criteria are set in accordance with the causal significance assessment of power transport. The array size of the distribution of causal significance indicators for chemical processes is also Nz×Ny×Nx.

[0106] Step 144: Perform spatial neighborhood smoothing filtering on the distribution of causal significance indicators of power transmission in the spatial grid dimension to reduce the impact of high-frequency random fluctuations in space on the indicator distribution, and obtain a smoothed causal significance distribution of power transmission. Perform spatial neighborhood smoothing filtering on the distribution of causal significance indicators of chemical processes in the spatial grid dimension to reduce the impact of high-frequency random fluctuations in space on the indicator distribution, and obtain a smoothed causal significance distribution of chemical processes.

[0107] Spatial neighborhood smoothing filtering employs a three-dimensional Gaussian filter kernel. The radii of the filter kernel in the vertical, latitude, and longitude dimensions are set to Rz, Ry, and Rx grid spacings, respectively, and the standard deviations of the Gaussian kernel in the three dimensions are set to σz, σy, and σx, respectively. For each grid location (z, y, x), the causal significance indices of all neighboring grids (z+dz, y+dy, x+dx) within its neighborhood that satisfy the conditions that the absolute value of the vertical offset is less than or equal to Rz, the absolute value of the latitude offset is less than or equal to Ry, and the absolute value of the longitude offset is less than or equal to Rx are weighted and averaged according to Gaussian weighting coefficients.

[0108] The Gaussian weighting coefficients w(dz, dy, dx) are calculated as follows: First, calculate the vertical weighting factor wz = exp(-dz) / dx. 2 ÷(2×σz 2 ), latitude direction weighting factor wy=exp(-dy) 2 ÷(2×σy 2 Longitude direction weighting factor wx = exp(-dx) 2 ÷(2×σx 2 Then, the three weighting factors are multiplied to obtain the final weighting coefficient w(dz, dy, dx) = wz × wy × wx. The weighted average calculation process is as follows: multiply the causal significance index of each grid in the neighborhood by its corresponding Gaussian weight coefficient, sum them, and then divide by the sum of all weight coefficients. The filtering process at the boundary grid position uses a mirror expansion method to fill the missing neighbor grid values, that is, when z + dz exceeds the effective range of the vertical layer, the filling value is obtained by reflecting the boundary layer as the axis of symmetry. The smoothed causal significance distribution after filtering has the same array size Nz × Ny × Nx as the original index distribution.

[0109] Step 145: Generate a set of power causal intensity distributions based on the mapping relationship between the smooth power transmission causal significance distribution and the preset index. The intensity value of each spatial grid in the power causal intensity distribution set corresponds to the significance level of the corresponding position in the smooth power transmission causal significance distribution in a monotonically increasing manner.

[0110] The preset index mapping relationship is a piecewise linear mapping function. Let the smoothed significance index value be S, and the mapped dynamic causal intensity value be I. When S is less than S1, I = k1 × S, where k1 is the slope of the first segment of the mapping. When S is between S1 and S2, I = I1 + k2 × (S - S1), where I1 is the mapping value when S equals S1, and k2 is the slope of the second segment of the mapping. When S is greater than S2, I = 1.0. The values ​​of the slopes k1 and k2 ensure that the mapping function is continuous at the segment points S1 and S2. This mapping relationship guarantees that the dynamic causal intensity value changes continuously and monotonically in the interval between 0 and 1. The above mapping calculation is performed on all Nz × Ny × Nx spatial grid locations to obtain the dynamic causal intensity distribution set, with an array size of Nz × Ny × Nx.

[0111] Step 146: Generate a set of chemical causal intensity distributions based on the mapping relationship between the causal significance distribution of smoothed chemical processes and preset indices. The intensity value of each spatial grid in the chemical causal intensity distribution set corresponds to the significance level of the corresponding position in the causal significance distribution of smoothed chemical processes in a monotonically increasing manner.

[0112] The chemical causal intensity distribution set is generated using the same piecewise linear mapping function as in step 145, with consistent parameter values ​​for the segmentation points S1 and S2 and the mapping slopes k1 and k2. The mapped chemical causal intensity values ​​also lie within the interval of 0 to 1, and the array size is Nz×Ny×Nx.

[0113] Step 150: Generate relative dominant role discrimination information based on the dynamic causal intensity distribution set and the chemical causal intensity distribution set.

[0114] Step 151: Perform a grid-by-grid numerical comparison between the dynamic causal intensity value at each spatial grid location in the dynamic causal intensity distribution set and the chemical causal intensity value at the same spatial grid location in the chemical causal intensity distribution set to obtain the dominant discriminant label for each spatial grid location.

[0115] For each spatial grid location (z, y, x), the intensity value Iadv is obtained from the dynamic causal intensity distribution set, and the intensity value Ichem is obtained from the chemical causal intensity distribution set. The difference between the two is calculated as Δ = Iadv - Ichem. The specific rules for grid-by-grid numerical comparison are as follows: when Δ is greater than Δhigh, the dominant discriminant label of the corresponding spatial grid location is determined to be the dynamic transport dominant type. When Δ is less than -Δhigh, the dominant discriminant label of the corresponding spatial grid location is determined to be the chemical process dominant type. When the absolute value of Δ is less than or equal to Δhigh, the dominant discriminant label of the corresponding spatial grid location is determined to be the equilibrium contribution type. Wherein Δhigh is the preset upper threshold of the equilibrium determination range. The dominant discriminant labels of all Nz×Ny×Nx spatial grids constitute the dominant discriminant label distribution.

[0116] Step 152: Obtain the seasonal oscillation phase period division information corresponding to the target atmospheric region. The seasonal oscillation phase period division information includes the convective active phase period marking, the convective suppression phase period marking, and the phase transition transition period marking.

[0117] The seasonal oscillation phase segmentation information is determined by bandpass filtering the average outward longwave radiation anomaly time series of the target atmospheric region, based on the sign of the filtering result and the zero-crossing position. First, the average outward longwave radiation anomaly area within the spatial range of the target atmospheric region is calculated, resulting in a one-dimensional time series of length T. Bandpass filtering is then applied to this time series, with the passband set to the frequency range corresponding to a 30-day to 60-day period, to extract the signal components at the intraseasonal oscillation scale.

[0118] In the filtered time series, consecutive periods with negative values ​​exceeding -1 standard deviation are marked as active convection phases. Consecutive periods with positive values ​​exceeding +1 standard deviation are marked as suppressed convection phases. Periods with filtered values ​​between -1 and +1 standard deviation, as well as transition periods between different phases, are marked as phase transition periods. Each time step is assigned one of these three phase phase labels, generating a phase phase label sequence of length T.

[0119] Step 153: Perform type and quantity statistics processing on the dominant discrimination labels corresponding to all time points marked as active convection phases, and obtain the proportion of the number of first spatial grids for the dominant power transport type, the proportion of the number of second spatial grids for the dominant chemical process type, and the proportion of the number of third spatial grids for the equilibrium contribution type under the active phase.

[0120] The statistical processing iterates through all time steps marked as active convection phases, and for each time step, the dominant discrimination label distribution corresponding to that time step is obtained. The number of grids of the power transport dominant type (Nadv_active(t), the number of grids of the chemical process dominant type (Nchem_active(t), and the number of grids of the balanced contribution type (Nequ_active(t))) in the dominant discrimination label distribution of that time step are counted respectively.

[0121] The number of grids for each type of time step in all active convection phases is summed to obtain the total number of grids for the dominant dynamic transport type (Sum_Nadv_active), the total number of grids for the dominant chemical process type (Sum_Nchem_active), and the total number of grids for the equilibrium contribution type (Sum_Nequ_active). The total effective number of grids in the active phase (Total_active) is equal to the sum of these three summed values. The proportion of grids in the first space is equal to Sum_Nadv_active ÷ Total_active, the proportion of grids in the second space is equal to Sum_Nchem_active ÷ Total_active, and the proportion of grids in the third space is equal to Sum_Nequ_active ÷ Total_active.

[0122] Step 154: Perform type and quantity statistics processing on the dominant discrimination labels corresponding to all time points belonging to the convection suppression phase period marker, and obtain the proportion of the number of fourth spatial grids of the dominant power transport type, the proportion of the number of fifth spatial grids of the dominant chemical process type, and the proportion of the number of sixth spatial grids of the equilibrium contribution type under the suppression phase.

[0123] The statistical processing method is the same as in step 153, but the statistical object is replaced with the time step corresponding to the convection suppression phase period. The number of grids of the dominant power transport type (Sum_Nadv_supp), the number of grids of the dominant chemical process type (Sum_Nchem_supp), and the number of grids of the equilibrium contribution type (Sum_Nequ_supp) are accumulated within the suppression phase period. The total effective number of grids in the suppression phase (Total_supp) is equal to the sum of the above three accumulated values. The proportion of the number of grids in the fourth space is equal to Sum_Nadv_supp ÷ Total_supp, the proportion of the number of grids in the fifth space is equal to Sum_Nchem_supp ÷ Total_supp, and the proportion of the number of grids in the sixth space is equal to Sum_Nequ_supp ÷ Total_supp.

[0124] Step 155: Perform type and quantity statistics processing on the dominant discrimination labels corresponding to all time points marked as phase transition periods, and obtain the proportion of the number of seventh spatial grids for the dominant power transport type, the proportion of the number of eighth spatial grids for the dominant chemical process type, and the proportion of the number of ninth spatial grids for the equilibrium contribution type under the phase transition.

[0125] The statistical processing method is the same as in step 153, but the statistical object is replaced with the time step corresponding to the phase transition period marker. The number of grids of the dominant power transport type (Sum_Nadv_tran), the number of grids of the dominant chemical process type (Sum_Nchem_tran), and the number of grids of the equilibrium contribution type (Sum_Nequ_tran) are accumulated separately during the phase transition period. The total effective number of grids in the phase transition (Total_tran) is equal to the sum of the above three accumulated values. The proportion of grids in the seventh space is equal to Sum_Nadv_tran ÷ Total_tran, the proportion of grids in the eighth space is equal to Sum_Nchem_tran ÷ Total_tran, and the proportion of grids in the ninth space is equal to Sum_Nequ_tran ÷ Total_tran.

[0126] Step 156: Generate relative dominant effect discrimination information based on the proportions of the number of grids in the first space, the second space, the fourth space, the fifth space, the seventh space, and the eighth space.

[0127] The information on relative dominance is organized in a structured text format, containing: the proportion of dynamic transport dominance in the active phase relative to the number of spatial grids in the first spatial grid, and the proportion of chemical processes dominance relative to the number of spatial grids in the second spatial grid; the proportion of dynamic transport dominance in the suppressive phase relative to the number of spatial grids in the fourth spatial grid, and the proportion of chemical processes dominance relative to the number of spatial grids in the fifth spatial grid; and the proportion of dynamic transport dominance in the transition phase relative to the number of spatial grids in the seventh spatial grid, and the proportion of chemical processes dominance relative to the number of spatial grids in the eighth spatial grid. This information quantitatively reveals the evolution of the relative importance of dynamic transport and chemical processes in regulating ozone distribution during different phases of the intraseasonal oscillation.

[0128] Following step 150, the method further includes: Step 210: Perform spatial grid mapping association processing on the dominant discrimination label corresponding to each spatial grid position in the relative dominant action discrimination information and the flow field direction distribution data in the convection parameter field set to generate a coupling representation unit between the dominant type and streamline direction at each spatial grid position, and construct a spatial coupling representation set based on the coupling representation units corresponding to all spatial grids.

[0129] The flow field direction distribution data is calculated using a streamline tracing algorithm on the horizontal wind field. For each spatial grid location and each time step, the zonal wind component U and the meridional wind component V are obtained for that grid location. Starting from that grid location, the streamline path is traced forward, with a tracing step size of Δs and a tracing number of steps of Nsteps. The formula for calculating the coordinates of the next location at each tracing step is: the new longitude equals the current longitude plus the U component divided by the horizontal wind vector magnitude multiplied by Δs divided by the Earth's radius divided by the cosine of the current latitude; the new latitude equals the current latitude plus the V component divided by the horizontal wind vector magnitude multiplied by Δs divided by the Earth's radius. The horizontal wind vector magnitude is equal to the square root of the sum of U² and V². After continuously tracing for Nstep steps, the average value of the direction vectors at each point on the path is taken as the representative value of the flow field direction at that grid location. The representative value of the flow field direction is a two-dimensional vector (Umean, Vmean). The representative flow field direction values ​​of all time steps are averaged over the time dimension to obtain the representative flow field direction data for that grid location.

[0130] The coupling representation unit is constructed as follows: the dominant discrimination label is converted into a 3D vector using one-hot encoding. The dominant type of power transport is encoded as (1, 0, 0), the dominant type of chemical process is encoded as (0, 1, 0), and the equilibrium contribution type is encoded as (0, 0, 1). This 3D one-hot vector is then concatenated with the two horizontal components (Umean, Vmean) of the flow field direction representative value to form a 5D coupling representation vector. The coupling representation vectors of all Nz×Ny×Nx grids are arranged spatially to form a spatial coupling representation set with an array size of Nz×Ny×Nx×5.

[0131] Step 220: Invoke the preset spatiotemporal backtracking processing mechanism to perform streamline reverse tracing processing from the current time segment to the historical time segment along the flow field direction on the spatial coupling representation set. Project the dominant discriminant label of each spatial grid position in the reverse direction along the corresponding streamline to the upstream spatial grid position corresponding to the time backtracking threshold to obtain the spatial position migration trajectory of each dominant discriminant label in the time backtracking direction.

[0132] The pre-defined spatiotemporal backtracking mechanism is implemented based on a semi-Lagrange backward trajectory tracking algorithm. The time backtracking threshold is set to Tback hours, and the corresponding time step is Nback = Tback ÷ Δt. For each spatial grid position in the current time segment t, the representative values ​​of the flow field direction (Umean, Vmean) at that position are obtained from the spatial coupling representation set. The direction of streamline reverse tracking is taken as (-Umean, -Vmean), indicating backtracking along the streamline in the reverse direction.

[0133] The reverse tracking calculation process is as follows: Starting from the current grid position's 3D coordinates (z, y, x), iterate Nback tracking steps. Each tracking step is calculated as follows: The reverse displacement vector is determined based on the zonal wind component U and the meridional wind component V at the current position. The new longitude coordinates are equal to the current longitude minus the U component divided by the horizontal wind vector magnitude multiplied by the tracking step size Δs divided by the Earth's radius divided by the cosine of the current latitude; the new latitude coordinates are equal to the current latitude minus the V component divided by the horizontal wind vector magnitude multiplied by the tracking step size Δs divided by the Earth's radius; the new vertical coordinates are determined based on the isentropic surface assumption. After each tracking step, the current position is updated to the new coordinates, and the spatial coordinates of that step are recorded. After Nback iterations, the upstream spatial grid position (z_back, y_back, x_back) is reached. The dominant discrimination label of this starting grid position is assigned to the tracking endpoint position. The spatial coordinates of each step in the tracking process are continuously recorded to form a spatial position sequence of length Nback+1. This sequence is the spatial position migration trajectory of the dominant discriminant label of the starting grid in the time backtracking direction.

[0134] Step 230: Perform spatial range extension analysis on all spatial location migration trajectories corresponding to the dominant power transmission type to extract the geometric extension span characteristics of the dominant power transmission path, and perform spatial range extension analysis on all spatial location migration trajectories corresponding to the dominant chemical process type to extract the geometric extension span characteristics of the dominant chemical process path. Generate dominant action propagation path description information based on the geometric extension span characteristics of the dominant power transmission path and the dominant chemical process path.

[0135] Spatial extent extension analysis is achieved by calculating the great circle distance between the starting and ending points of each migration trajectory. For a migration trajectory with starting coordinates (lat1, lon1) and ending coordinates (lat2, lon2), the calculation logic of the great circle distance D is as follows: First, calculate the latitude difference Δlat and longitude difference Δlon between the two points. Then, calculate the intermediate variable a, which equals sin²(Δlat÷2) plus cos(lat1)×cos(lat2)×sin²(Δlon÷2). Next, calculate the intermediate variable c, which equals 2×atan²(√a, √(1-a)). Finally, the great circle distance D equals the Earth's radius multiplied by c.

[0136] For all spatial grids with the dominant identification label "Power Transmission Dominant Type," obtain the spatial location migration trajectory for each path and calculate its great circle distance. Sort the great circle distance values ​​of all power transmission dominant type trajectories in ascending order, and take the distance value corresponding to the 50th percentile after sorting as the typical extension span Dadv_50 of the power transmission dominant path. Take the distance value corresponding to the 90th percentile as the maximum extension span Dadv_90 of the power transmission dominant path.

[0137] The same operation was performed on all spatial grids labeled as chemical process-dominant to obtain the typical span (Dchem_50) and maximum span (Dchem_90) of the chemical process-dominant path. The propagation path description information of the dominant effect includes the above four span values, characterizing the differences in the spatial influence range between dynamic transport-dominant and chemical process-dominant effects.

[0138] Following step 150, the method further includes: Step 310: Obtain the seasonal oscillation phase attribute label sequence of the target atmospheric region in multiple consecutive time segments. Align the dominant discrimination label of each spatial grid position in the relative dominant action discrimination information with the seasonal oscillation phase attribute label sequence in different time segments and perform time axis alignment and spatial grid correspondence processing to obtain the dominant discrimination label time series of each spatial grid position under the phase attribute label constraint.

[0139] Multiple consecutive time segments are divided into M time segments, the length of which can be set according to the analysis requirements. For each time segment, the seasonal oscillation phase attribute markers corresponding to that time segment are obtained according to the method in step 152, generating a phase attribute marker sequence of length M. The time axis alignment process matches the phase attribute marker sequence with the dominant discriminant label of each spatial grid position according to the timestamp information of each time segment.

[0140] After matching, each spatial grid location yields a time series of length M. Each element in the series is a binary tuple. The first component of the binary tuple is the intraseasonal oscillation phase attribute label for that time segment, and the second component is the dominant discriminant label type for that time segment. The binary tuple data structure contains two fields: the label phase type is active phase, suppressed phase, or transition phase; the dominant label type is dominated by power transport, chemical process, or equilibrium contribution.

[0141] Step 320: Perform cross-phase type transfer statistical processing on the dominant discrimination label time series corresponding to each spatial grid location to generate first transfer description information characterizing the frequency of the transition between different phase attributes of the dominant state of power transport and second transfer description information characterizing the frequency of the transition between different phase attributes of the dominant state of chemical process. Combine the first transfer description information and the second transfer description information into a phase transfer association description set.

[0142] Cross-phase type transition statistics are achieved by constructing a state transition matrix. The row index of the state transition matrix is ​​the phase attribute of the starting time segment, and the column index is the phase attribute of the ending time segment. The phase attributes include three types: active phase, suppressed phase, and transition phase. Therefore, the transition matrix is ​​a 3x3 square matrix.

[0143] For the dominant state of power transmission, iterate through two consecutive time segments in the dominant discrimination label time series. If the dominant discrimination label of the previous time segment is power transmission dominant type, then obtain the phase attribute of the previous time segment as the starting state and the phase attribute of the next time segment as the ending state, and increment 1 at the corresponding row and column positions of the transition matrix. After traversing all consecutive time segment pairs, the frequency transition matrix of the dominant state of power transmission is obtained. Normalize each element in the matrix by dividing it by the sum of all elements in the matrix to obtain the first transition description information, which is a 3×3 probability transition matrix. The element in the i-th row and j-th column of the matrix represents the probability of transitioning from the i-th phase to the j-th phase under the dominant state of power transmission.

[0144] Perform the same operation on the dominant state of the chemical process: traverse two consecutive time segments in the dominant discrimination label time series. If the dominant discrimination label of the previous time segment is the dominant type of the chemical process, increment the corresponding phase attribute transition position in the transition matrix by 1. After normalization, the second transition description information is obtained, which is also a 3×3 probability transition matrix.

[0145] The first and second transition description information are combined into a phase transition association description set, which contains two 3×3 probability transition matrices.

[0146] Step 330: Perform spatial clustering analysis on the transfer description information corresponding to each spatial grid position in the phase transfer association description set, and generate phase-dependent ozone control mechanism description information based on the distribution characteristics of the clustering regions of transfer frequency in the spatial dimension.

[0147] Spatial clustering analysis employs local spatial autocorrelation statistics. For each phase transition type in the transition description information, the transition probability value corresponding to that type of transition at each spatial grid location is extracted from the phase transition association description set, forming a spatial transition probability distribution field with an array size of Nz×Ny×Nx.

[0148] For each spatial grid location, the local Moran index is calculated. The local Moran index I is calculated as follows: First, the deviation of the transition probability value at that grid location from the average transition probability value for the entire region is calculated. Then, for each neighboring grid within the grid's neighborhood, the deviation of the neighboring grid's transition probability value from the average value for the entire region is calculated. These deviations are multiplied, multiplied by a spatial weighting coefficient, summed over all neighboring grids, and then multiplied by a normalization factor. The spatial weighting coefficient uses distance decay weighting, with closer neighboring grids having a larger weight. The normalization factor is the variance of the entire region's transition probability.

[0149] The significance level of the local Moran's index was determined using a 999-permutation test. The process of the random permutation test is as follows: keeping the spatial grid positions unchanged, the transition probability values ​​of each grid are randomly shuffled 999 times. Each time, the local Moran's index for each grid is calculated, resulting in 999 randomly distributed samples. The calculated local Moran's index is compared with the randomly distributed samples, and the percentile rank of the actual value in the random distribution is calculated. If the actual value is greater than 95% of the random sample values, the clustering characteristic of that grid position is considered statistically significant, with a Z-value greater than 1.96.

[0150] Grid locations with significantly positive local Moran's indices and Z-values ​​greater than 1.96 are identified as high-value clustering regions, indicating a spatially high-value clustering distribution of this type of phase transition probability in these areas. Grid locations with significantly negative local Moran's indices and Z-values ​​less than -1.96 are identified as low-value anomaly regions. The phase-dependent ozone control mechanism description information correlates the spatial distribution characteristics of the above-mentioned clustering regions with the corresponding phase transition paths, clarifying the spatial distribution patterns of dynamic transport-dominated or chemical process-dominated phase transition stages. For example, it identifies the latitude and longitude ranges and vertical height levels where high-probability clustering regions dominated by dynamic transport are mainly distributed during the transition from an active phase to an inhibited phase.

[0151] This application embeds the constraints of the Navier-Tox equation and the linear theory of gravity waves into a pre-defined physical information network as physical regularization constraints. During the spatiotemporal evolution physical constraint processing of the convective parameter field set and the gravity wave activity field set, it achieves the decoupling and separation of the power transport process and the chemical process under the guidance of physical laws. This ensures that the output decoupled power transport contribution field set and decoupled chemical process contribution field set possess physical consistency and interpretability. Furthermore, by performing time-series correlation analysis on the decoupled contribution field set and the ozone concentration field set respectively, it extracts the time-series response feature sets between the power transport contribution time series and the chemical process contribution time series and the ozone concentration time series, effectively capturing the leading and lagging response relationships and attenuation characteristics of different physical processes to ozone concentration changes in the time dimension. Next, the two types of time-series response feature sets are input into a pre-defined causal inference model for causal intensity inference processing, obtaining the power causal intensity distribution set and the chemical causal intensity distribution set respectively, thereby quantitatively characterizing the spatial distribution of the causal influence intensity of power transport and chemical processes on ozone concentration changes. Finally, based on the two sets of causal intensity distribution sets, relative dominant role discrimination information was generated, which enabled an objective determination of the relative importance of the contribution of power transport and chemical process at different spatial locations, and improved the accuracy and reliability of ozone budget attribution analysis near the tropopause.

[0152] It should be noted that although the decomposition logic of the dynamic transport contribution in the technical solution of this application involves the joint decoupling of convective effective potential energy and gravitational fluctuation flux, its mathematical implementation falls within the scope of numerical solutions that can be conventionally configured based on atmospheric dynamic equations by those skilled in the art. Regarding the separation of the advection transport term and the divergence term in the ozone continuity equation, those skilled in the art know that in atmospheric chemical transport models, a divergence operator in flux form is typically used for characterization; that is, the dynamic tendency term can be expressed as the divergence of the negative product of ozone concentration and wind field vector. When faced with the need to physically decompose a single scalar tendency value output by a neural network, the Helmholtz decomposition theorem in computational fluid dynamics, or the use of the horizontal wind field divergence field and advection gradient field pre-calculated based on reanalysis data as spatial basis functions, can be used to solve the combination coefficients through local linear regression to approximately characterize the relative contribution weights of flux divergence and advection transport at that grid point. This is a conventional fitting method in model post-processing and feature attribution analysis, and its accuracy is limited by the input field resolution but does not present any fundamental implementation obstacles.

[0153] Regarding the issue of momentum constraint weight decay caused by the Coriolis parameter approaching zero in the low-latitude equatorial region, those skilled in the art are well aware that, within the framework of tropical atmospheric dynamics, the applicable form of the Navier-Stokes equations is usually transformed into shallow-water equations or constraints are constructed using scale analysis based on Rossby number or Richardson number. In the construction of regularization terms for actual network training, the constraint form can be dynamically adjusted according to the latitudinal zone of the target region. For example, in the extratropical region, the focus is on geostrophic balance and momentum flux divergence constraints, while in the equatorial waveguide region, it automatically degenerates into a physical regularization strategy dominated by the continuity equation and gravity wave dispersion relation. Alternatively, a specific momentum constraint expression under the equatorial beta plane approximation can be introduced to replace the Coriolis parameter weight correction. Such approaches have mature application precedents in large-scale circulation diagnosis and parameterization scheme evaluation.

[0154] Regarding the consistency of wavenumber definition and dimensions in gravity wave dispersion relations, those skilled in the art, with general knowledge of atmospheric wave physics, can clearly distinguish the conversion relationship between angular wavenumber and cyclic wavenumber in the expression. When referring to the calculation of the vertical wavenumber square, it actually refers to the process of converting the wavenumber components corresponding to the spectral peaks obtained from the Fast Fourier Transform according to the Nyquist sampling theorem and the vertical grid spacing into physical wavenumbers. In standard atmospheric science calculation conventions, the angular wavenumber vector with the correct length dimension reciprocal is obtained by multiplying the dimensionless wavenumber index by twice pi and dividing by the physical length of the vertical domain. This vector is then substituted into the theoretical dispersion polynomial of the natural frequency, buoyancy frequency, and wavenumber vector for deviation evaluation. This type of wavenumber spectrum analysis and dispersion feature extraction process has standardized implementation references in gravity wave parameter inversion from observational data and gravity wave source-sink diagnosis in high-resolution numerical models, ensuring logical consistency of physical constraints at the dimensional matching and numerical calculation levels.

[0155] Based on the same inventive concept, embodiments of this application also provide a convective gravity wave decoupling server. See also... Figure 2 As shown, this is a schematic diagram of a possible convective-gravity wave decoupling server provided in an embodiment of this application. Figure 2 In the process, the convective gravity wave decoupling server 200 includes a processor 210 and a memory 220. The processor 210 and the memory 220 are interconnected via a communication bus. The memory 220 stores computer programs that can be executed by the processor 210. By executing the instructions stored in the memory 220, the processor 210 can perform the steps of the aforementioned convective gravity wave decoupling method based on a physical information neural network.

[0156] Based on the same inventive concept, embodiments of this application provide a computer-readable storage medium including a computer program. When the computer program runs on a convective gravity wave decoupling server, it causes the server to execute the steps of the aforementioned convective gravity wave decoupling method based on a physical information neural network. In some possible implementations, various aspects of the convective gravity wave decoupling method based on a physical information neural network provided in this application can also be implemented as a program product, including a computer program. When the program product runs on a convective gravity wave decoupling server, it causes the server to execute the steps of the aforementioned convective gravity wave decoupling method based on a physical information neural network. For example, the convective gravity wave decoupling server can execute... Figure 1 The steps are shown in the diagram. The computer-readable storage medium includes volatile or non-volatile or a combination thereof, and may be removable or non-removable. Examples of computer-readable storage media include, but are not limited to, phase-change random access memory (PRAM), static random access memory (SRAM), dynamic random access memory (DRAM), other types of random-access memory (RAM), read-only memory (ROM), electrically erasable programmable read-only memory (EEPROM), erasable programmable read-only memory (EPROM), programmable read-only memory (PROM), flash memory or other memory technologies, CD-ROM, Digital Video Disc (DVD) or other optical storage, magnetic tape, magnetic disk storage or other magnetic storage devices, or any other non-transfer medium.

[0157] like Figure 3The diagram shown is a functional block diagram of the convective gravity wave decoupling server provided in this embodiment of the application. The convective gravity wave decoupling server includes a convective gravity wave decoupling device, which includes: The data field acquisition module is used to acquire the set of convective parameter fields, gravity wave activity fields, and ozone concentration fields of the target atmospheric region. The physical constraint processing module is used to input the set of convection parameter fields and the set of gravity wave activity fields into a preset physical information network for spatiotemporal evolution physical constraint processing, and output the set of decoupled power transport contribution fields and the set of decoupled chemical process contribution fields. The spatiotemporal evolution physical constraint processing uses the Navier-Tex equation constraint terms and the gravity wave linear theory constraint terms as physical regularization constraints in the network training process. The time-series correlation analysis module is used to perform a first time-series correlation analysis on the decoupled power transmission contribution field set and the ozone concentration field set to extract a first time-series response feature set between the power transmission contribution time series and the ozone concentration time series, and to perform a second time-series correlation analysis on the decoupled chemical process contribution field set and the ozone concentration field set to extract a second time-series response feature set between the chemical process contribution time series and the ozone concentration time series; The causal intensity inference module is used to input the first time-series response feature set and the second time-series response feature set into a preset causal inference model to perform dynamic causal intensity inference processing to obtain a dynamic causal intensity distribution set, and to perform chemical causal intensity inference processing on the second time-series response feature set through the preset causal inference model to obtain a chemical causal intensity distribution set; The discrimination information generation module is used to generate relative dominant effect discrimination information based on the dynamic causal intensity distribution set and the chemical causal intensity distribution set.

[0158] Accordingly, this application also provides a computer program product, which includes a computer program or instructions that, when executed by a processor, cause the processor to implement the steps in the above method embodiments. It should be understood that each step or combination of steps in the above method flow can be implemented by the computer program or instructions. Furthermore, these computer programs or instructions can be applied to the processor of a general-purpose computer, a special-purpose computer, an embedded processor, or other programmable data processing device, enabling the processor of the general-purpose computer, special-purpose computer, embedded processor, or other programmable data processing device to function as an apparatus for implementing the corresponding functions in the above method embodiments.

[0159] Those skilled in the art will clearly understand that, for the sake of convenience and brevity, the specific working processes of the systems, devices, and units described above can be referred to the corresponding processes in the foregoing method embodiments, and will not be repeated here.

[0160] It should also be noted that the terms "comprising," "including," or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or apparatus that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such process, method, article, or apparatus. Unless otherwise specified, an element defined by the phrase "comprising one..." does not exclude the presence of other identical elements in the process, method, article, or apparatus that includes that element.

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

Claims

1. A method for decoupling convective-gravity waves based on a physical information neural network, characterized in that, include: Acquire the set of convective parameter fields, gravity wave activity fields, and ozone concentration fields for the target atmospheric region; The set of convection parameters and the set of gravity wave activity fields are input into a preset physical information network for spatiotemporal evolution physical constraint processing, and the decoupled power transport contribution field set and the decoupled chemical process contribution field set are output. The spatiotemporal evolution physical constraint processing uses the Navier-Tex equation constraint term and the gravity wave linear theory constraint term as the physical regularization restriction condition in the network training process. The decoupled power transmission contribution field set and the ozone concentration field set are subjected to a first time series correlation analysis to extract a first time series response feature set between the power transmission contribution time series and the ozone concentration time series. The decoupled chemical process contribution field set and the ozone concentration field set are subjected to a second time series correlation analysis to extract a second time series response feature set between the chemical process contribution time series and the ozone concentration time series. The first time-series response feature set and the second time-series response feature set are input into a preset causal inference model to perform dynamic causal intensity inference processing to obtain a dynamic causal intensity distribution set, and the second time-series response feature set is then processed by the preset causal inference model to perform chemical causal intensity inference processing to obtain a chemical causal intensity distribution set; Relative dominance discrimination information is generated based on the set of dynamic causal intensity distributions and the set of chemical causal intensity distributions.

2. The method according to claim 1, characterized in that, The step of inputting the set of convective parameter fields and the set of gravity wave activity fields into a preset physical information network for spatiotemporal evolution physical constraint processing, and outputting a decoupled power transport contribution field set and a decoupled chemical process contribution field set, includes: Extract convective effective potential energy distribution units at multiple spatial grid locations from the convective parameter field set, and extract gravitational wave flux units at multiple spatial grid locations from the gravitational wave activity field set; The convective effective potential energy distribution unit and the gravitational fluctuation flux unit are aligned with spatial grid coordinates to obtain aligned convective effective potential energy distribution and aligned gravitational fluctuation flux distribution. The aligned convective effective potential energy distribution and the aligned gravitational fluctuation flux distribution are sequentially combined along the time axis to generate a joint spatiotemporal grid data set. The shared feature extraction layer in the preset physical information network is invoked to perform spatial dimension downsampling and temporal dimension convolution processing on the joint spatiotemporal grid data set to obtain a spatiotemporal coupled potential feature representation set. The power transmission mapping branch in the preset physical information network is invoked to perform grid-by-grid feature transformation processing on the spatiotemporal coupled potential feature representation set to generate a preliminary distribution set of power transmission contributions. The chemical process mapping branch in the preset physical information network is invoked to perform grid-by-grid feature transformation processing on the spatiotemporal coupled potential feature representation set to generate a preliminary distribution set of chemical process contributions. Based on the preliminary distribution set of the power transmission contribution and the preliminary distribution set of the chemical process contribution, a physical rule compliance evaluation function is constructed. The physical rule compliance evaluation function includes a compliance evaluation sub-function of the Navier-Tox equation and a compliance evaluation sub-function of the linear theory of gravity waves. A joint optimization objective function is constructed based on the output value of the physical rule compliance evaluation function and the data reconstruction difference value. The internal connection weights of the preset physical information network are iteratively adjusted according to the joint optimization objective function. When the iteration adjustment process reaches a preset stable state, the preliminary distribution set of power transmission contribution output by the power transmission mapping branch is output as the decoupled power transmission contribution field set, and the preliminary distribution set of chemical process contribution output by the chemical process mapping branch is output as the decoupled chemical process contribution field set.

3. The method according to claim 2, characterized in that, The construction of a physical rule compliance evaluation function based on the preliminary distribution set of power transmission contributions and the preliminary distribution set of chemical process contributions includes: From the preliminary distribution set of power transport contribution, the velocity field divergence component distribution unit and the momentum field advection component distribution unit are separated; Based on the analysis of the divergence component distribution unit of the velocity field and the divergence deviation distribution of the incompressible fluid continuity condition, the continuity constraint compliance component of the Navier-Tox equation is obtained. Based on the analysis of the momentum field advection component distribution unit and the momentum balance deviation distribution between the momentum field and the gravity fluctuation source and sink terms, the momentum constraint compliance component of the Navier-Tox equation is obtained. The vertical wavenumber characteristic distribution and the horizontal wavenumber characteristic distribution of gravity waves are separated from the gravity fluctuation flux unit of the gravity wave activity field set, and the background flow field adjustment component distribution caused by the power transmission process is extracted from the preliminary distribution set of power transmission contribution. Based on the vertical wavenumber characteristic distribution of gravity waves, the horizontal wavenumber characteristic distribution of gravity waves, and the background flow field adjustment component distribution, the dispersion relation response deviation distribution and the refraction deviation distribution of energy propagation direction of gravity waves are analyzed, and the dispersion constraint compliance component and the refraction constraint compliance component of the linear theory of gravity waves are obtained. The continuity constraint compliance component and the momentum constraint compliance component are spatially averaged and fused to generate the compliance evaluation sub-function of the Navier-Tox equation. The dispersion constraint compliance component and the refractive constraint compliance component are spatially averaged and fused to generate a compliance evaluation sub-function for the linear theory of gravity waves. The conformity evaluation sub-function of the Navier-Tox equation and the conformity evaluation sub-function of the linear theory of gravity waves are weighted and summed to generate the output value of the conformity evaluation function of the physical rule.

4. The method according to claim 1, characterized in that, The step of performing a first time-series correlation analysis on the decoupled power transport contribution field set and the ozone concentration field set to extract a first time-series response feature set between the power transport contribution time series and the ozone concentration time series, and performing a second time-series correlation analysis on the decoupled chemical process contribution field set and the ozone concentration field set to extract a second time-series response feature set between the chemical process contribution time series and the ozone concentration time series, includes: For each spatial grid in the decoupled power transmission contribution field set, extract time series units of power transmission contribution intensity arranged along the time direction; for the same spatial grid in the ozone concentration field set, extract time series units of ozone concentration value arranged along the time direction; perform synchronous sliding time window division processing on the time series units of power transmission contribution intensity and the time series units of ozone concentration value to obtain multiple sets of time window corresponding to power transmission contribution window segments and ozone concentration window segments; For each time window, a lag cross-correlation analysis is performed on the power transport contribution window segment and the ozone concentration window segment to generate the first lag correlation coefficient change curve at each spatial grid location. The peak points of the correlation coefficient in the positive lag region and the negative lag region are identified in the first lag correlation coefficient change curve, and the first lead response time parameter corresponding to the positive lag peak point and the first lag response time parameter corresponding to the negative lag peak point are extracted. Analyze the correlation coefficient decay trend on both sides of the peak point of the first lag correlation coefficient change curve as the absolute value of the lag increases, and extract the first correlation coefficient decay characteristic parameter that characterizes the rate of decay. The first lead response time parameter, the first lag response time parameter and the first correlation coefficient decay characteristic parameter are combined to form the first time series response feature set corresponding to each spatial grid. Extract time series units of chemical process contribution intensity arranged along the time direction from each spatial grid in the decoupled chemical process contribution field set; perform synchronous sliding time window division on the time series units of chemical process contribution intensity and the time series units of ozone concentration value to obtain multiple sets of time window corresponding to chemical process contribution window segments and ozone concentration window segments. For each time window, a hysteresis cross-correlation analysis is performed on the chemical process contribution window segment and the ozone concentration window segment to generate a second hysteresis correlation coefficient change curve at each spatial grid location. The peak points of the correlation coefficient in the positive hysteresis region and the peak points of the correlation coefficient in the negative hysteresis region are identified in the second hysteresis correlation coefficient change curve. The second lead response time parameter corresponding to the positive hysteresis peak point and the second hysteresis response time parameter corresponding to the negative hysteresis peak point are extracted. Analyze the correlation coefficient decay trend on both sides of the peak point of the second lag correlation coefficient change curve as the absolute value of the lag increases, and extract the second correlation coefficient decay characteristic parameter that characterizes the rate of decay. The second lead response time parameter, the second lag response time parameter, and the second correlation coefficient decay characteristic parameter are combined to form the second time series response feature set corresponding to each spatial grid.

5. The method according to claim 1, characterized in that, The step of inputting the first time-series response feature set and the second time-series response feature set into a preset causal inference model, calling the preset causal inference model to perform dynamic causal intensity inference processing on the first time-series response feature set to obtain a dynamic causal intensity distribution set, and performing chemical causal intensity inference processing on the second time-series response feature set through the preset causal inference model to obtain a chemical causal intensity distribution set includes: The first response feature vector corresponding to each spatial grid in the first time-series response feature set is concatenated with the second response feature vector corresponding to the same spatial grid in the second time-series response feature set to generate a fusion causal feature representation set of multiple spatial grids; The first causal test branch in the preset causal inference model is invoked to perform a causal significance assessment of the fused causal feature representation set, analyze the degree of improvement in the ozone concentration time series prediction ability before and after the introduction of the power transmission contribution time series, and obtain the distribution of power transmission causal significance index. The second causal test branch in the preset causal inference model is invoked to perform chemical process causal significance assessment on the fused causal feature representation set. The degree of improvement in the prediction ability of ozone concentration time series before and after the introduction of chemical process contribution time series is analyzed, and the distribution of chemical process causal significance index is obtained. The causal significance index of power transmission is processed by spatial neighborhood smoothing filtering in the spatial grid dimension to reduce the influence of high-frequency random fluctuations in space on the index distribution, thus obtaining a smooth causal significance distribution of power transmission. The causal significance index of chemical process is processed by spatial neighborhood smoothing filtering in the spatial grid dimension to reduce the influence of high-frequency random fluctuations in space on the index distribution, thus obtaining a smooth causal significance distribution of chemical process. A set of power causal intensity distributions is generated based on the mapping relationship between the smooth power transmission causal significance distribution and the preset index. The intensity value of each spatial grid in the set of power causal intensity distributions corresponds to the significance level of the corresponding position in the smooth power transmission causal significance distribution in a monotonically increasing manner. A chemical causal intensity distribution set is generated based on the mapping relationship between the causal significance distribution of the smoothed chemical process and the preset index. The intensity value of each spatial grid in the chemical causal intensity distribution set corresponds to the significance level of the corresponding position in the causal significance distribution of the smoothed chemical process in a monotonically increasing manner.

6. The method according to claim 5, characterized in that, The first causal test branch in the preset causal inference model is invoked to perform a causal significance assessment of the fused causal feature representation set, analyzing the improvement in the predictive ability of ozone concentration time series before and after the introduction of the power transportation contribution time series, and obtaining the distribution of power transportation causal significance index, including: From the fused causal feature representation set, the historical observation sequence of power transmission contribution and the historical observation sequence of ozone concentration are decomposed. Based on the historical observation sequence of ozone concentration, a restricted prediction model that only relies on the historical information of ozone concentration itself is constructed. The current ozone concentration is forward extrapolated through the restricted prediction model to obtain the restricted ozone concentration extrapolation sequence. Obtain the degree of difference between the restricted ozone concentration projection sequence and the actual ozone concentration observation sequence; construct an unrestricted prediction model containing external input based on the historical observation sequence of power transport contribution and the historical observation sequence of ozone concentration; perform forward projection processing on the current ozone concentration through the unrestricted prediction model to obtain the unrestricted ozone concentration projection sequence. Obtain the unrestricted ozone concentration projection sequence and the actual ozone concentration observation sequence as an unrestricted projection difference index; construct a power transport causal significance evaluation function based on the relative change amplitude between the restricted projection difference index and the unrestricted projection difference index, and obtain the power transport causal significance evaluation value for each spatial grid. The causal significance assessment value of power transmission is compared with the preset significance discrimination conditions, and the distribution of causal significance index of power transmission is determined based on the comparison results.

7. The method according to claim 5 or 6, characterized in that, The second causal test branch in the preset causal inference model is invoked to perform a chemical process causal significance assessment on the fused causal feature representation set. The improvement in the predictive ability of ozone concentration time series before and after the introduction of the chemical process contribution time series is analyzed, resulting in the distribution of chemical process causal significance indicators, including: From the fused causal feature representation set, the historical observation sequence of chemical process contribution and the historical observation sequence of ozone concentration are decomposed; Based on the historical ozone concentration observation sequence, a restricted prediction model that relies solely on the historical information of the ozone concentration itself is constructed. The current ozone concentration is then forward-engineered using the restricted prediction model to obtain a restricted ozone concentration projection sequence. Obtain the index of the degree of difference between the restricted ozone concentration projection sequence and the actual ozone concentration observation sequence; construct an unrestricted prediction model containing external input based on the historical observation sequence of chemical process contribution and the historical observation sequence of ozone concentration; perform forward projection processing on the current ozone concentration through the unrestricted prediction model to obtain the unrestricted ozone concentration projection sequence. Obtain the unrestricted inference index between the unrestricted ozone concentration projection sequence and the actual ozone concentration observation sequence; construct a chemical process causal significance evaluation function based on the relative change amplitude between the restricted inference index and the unrestricted inference index, and obtain the chemical process causal significance evaluation value for each spatial grid. The causal significance assessment value of the chemical process is compared with the preset significance discrimination conditions, and the distribution of causal significance index of the chemical process is determined based on the comparison results.

8. The method according to any one of claims 1-6, characterized in that, The step of generating relative dominance discrimination information based on the dynamic causality intensity distribution set and the chemical causality intensity distribution set includes: The dynamic causal intensity value at each spatial grid position in the dynamic causal intensity distribution set is compared with the chemical causal intensity value at the same spatial grid position in the chemical causal intensity distribution set to obtain the dominant discriminant label for each spatial grid position. When the dynamic causal intensity value exceeds the chemical causal intensity value, the dominant discrimination label of the corresponding spatial grid position is determined as the dynamic transport dominant type; when the chemical causal intensity value exceeds the dynamic causal intensity value, the dominant discrimination label of the corresponding spatial grid position is determined as the chemical process dominant type; when the absolute difference between the dynamic causal intensity value and the chemical causal intensity value is less than the preset equilibrium judgment range, the dominant discrimination label of the corresponding spatial grid position is determined as the equilibrium contribution type. Acquire the intraseasonal oscillation phase period division information corresponding to the target atmospheric region. The intraseasonal oscillation phase period division information includes the convective active phase period marker, the convective suppression phase period marker, and the phase transition transition period marker. Perform type and quantity statistics processing on the dominant discrimination labels corresponding to all time points belonging to the active convection phase period marker to obtain the proportion of the first spatial grid number of the dominant power transport type, the proportion of the second spatial grid number of the dominant chemical process type, and the proportion of the third spatial grid number of the equilibrium contribution type under the active phase. Perform type and quantity statistics processing on the dominant discrimination labels corresponding to all time points belonging to the convection suppression phase period marker to obtain the proportion of the number of fourth spatial grids of the dominant power transport type, the proportion of the number of fifth spatial grids of the dominant chemical process type, and the proportion of the number of sixth spatial grids of the equilibrium contribution type under the suppression phase. Perform type and quantity statistics on the dominant discrimination labels corresponding to all time points belonging to the phase transition period marker to obtain the proportion of the number of seventh spatial grids of the dominant power transport type, the proportion of the number of eighth spatial grids of the dominant chemical process type, and the proportion of the number of ninth spatial grids of the equilibrium contribution type under the phase transition. Relative dominance discrimination information is generated based on the proportions of the first, second, fourth, fifth, seventh, and eighth spatial grids.

9. A convective-gravity wave decoupling server, characterized in that, It includes a processor and a memory, wherein the memory stores a computer program that, when executed by the processor, causes the processor to perform the steps of the convective gravity wave decoupling method based on a physical information neural network as described in any one of claims 1 to 8.

10. A computer-readable storage medium, characterized in that, The computer-readable storage medium includes a computer program that, when run on a convective gravity wave decoupling server, causes the convective gravity wave decoupling server to perform the steps of the convective gravity wave decoupling method based on a physical information neural network as described in any one of claims 1 to 8.

Images