A wide-area atmospheric pollutant concentration probabilistic prediction method, device and medium

Abstract

The application provides a wide-area atmospheric pollutant concentration probabilistic prediction method, equipment and medium, which is realized through a neural network model containing a three-dimensional space-time attention encoder, a mean prediction head and a covariance prediction head. The method first acquires a historical time input space-time cube containing multiple time steps, multiple spatial positions and multiple atmospheric parameters; inputs the space-time cube into the three-dimensional space-time attention encoder to extract a space-time feature map; then generates a future time deterministic prediction mean field through the mean prediction head and a prediction uncertainty covariance matrix through the covariance prediction head; and the two jointly define a multivariate Gaussian distribution of the future atmospheric state. The method realizes wide-area scale probabilistic prediction and uncertainty quantification, and provides accurate and reliable risk assessment basis for regional pollution prevention and control.

Description

Technical Field

[0001] This invention relates to the field of atmospheric environment prediction technology, and in particular to a method, equipment and medium for probabilistic forecasting of wide-area atmospheric pollutant concentrations. Background Technology

[0002] Atmospheric environment prediction is a key technology in the fields of environmental protection and public health. Traditional methods mainly rely on physical mechanism-based numerical simulations (such as WRF-Chem and CMAQ) to simulate the generation, diffusion, and transport processes of pollutants by solving complex systems of partial differential equations. Although these methods have a solid physical foundation, their computational cost is enormous, with a single high-resolution forecast often taking several hours, making it difficult to meet the needs of real-time emergency response. Furthermore, model optimization relies on manual parameter tuning, making it difficult to adaptively learn from massive amounts of historical observation data, leading to optimization bottlenecks. In recent years, data-driven deep learning methods (such as LSTM, Transformer, and CNN) have made breakthroughs in time series prediction and image processing, bringing new ideas to atmospheric environment prediction. However, existing deep learning methods still face many challenges: difficulties in fusing multi-source heterogeneous data, static spatial correlation modeling, tight coupling of time series prediction modules, and lack of uncertainty quantification, which limit their application in practical operations.

[0003] Uncertainty in wide-area forecasting stems from complex sources, including meteorological input errors, model approximation errors, and atmospheric chaotic characteristics. Existing mainstream methods (such as GraphCast and the Pangu meteorological model) primarily focus on deterministic forecasts, outputting single numerical values ​​and failing to quantify the reliability of prediction results, thus providing insufficient support for risk decision-making. Furthermore, prediction errors at different spatial locations often exhibit correlations (e.g., a region may have a systematically higher error rate), but current technologies struggle to model such spatially correlated uncertainties. Directly outputting the complete covariance matrix is ​​computationally too intensive and impractical for engineering implementation. Therefore, an efficient probabilistic forecasting method is needed that can output the uncertainty at each location, characterize spatial correlations, and ensure computational feasibility. Summary of the Invention

[0004] To address the aforementioned technical problems, the technical solution adopted by this invention is as follows:

[0005] According to a first aspect of the present invention, a probabilistic forecasting method for wide-area atmospheric pollutant concentrations is provided. The method is implemented through a pre-constructed neural network model, which includes a three-dimensional spatiotemporal attention encoder, a mean prediction head, and a covariance prediction head. The method includes the following steps:

[0006] The input spatiotemporal cube for obtaining historical time steps is obtained. The input spatiotemporal cube contains multiple time steps, multiple spatial locations, and multiple atmospheric parameters.

[0007] The input spatiotemporal cube is fed into an encoder based on a three-dimensional spatiotemporal attention mechanism to extract spatiotemporal feature maps.

[0008] Based on the spatiotemporal feature map, a deterministic forecast mean field for future times is generated using a mean prediction head.

[0009] Based on the spatiotemporal feature map, a covariance matrix for forecasting uncertainty is generated using a covariance prediction head.

[0010] The deterministic forecast mean field and the covariance matrix together define the multivariate Gaussian distribution of the atmospheric state at future times.

[0011] According to a second aspect of the present invention, an electronic device is provided, including a processor and a memory; the processor executes the steps of the method described in the first aspect of the present invention by invoking a program or instructions stored in the memory.

[0012] According to a third aspect of the present invention, a computer-readable storage medium is provided that stores a program or instructions that cause a computer to perform the steps of the method described in the first aspect of the present invention.

[0013] The present invention has at least the following beneficial effects:

[0014] This invention, building upon wide-area deterministic forecasting, adds a covariance prediction head and employs a covariance modeling method combining block diagonalization and low-rank global correlation to quantify the uncertainty in wide-area atmospheric environmental factor forecasts. This overcomes the shortcomings of existing wide-area forecast models that only output single numerical values ​​and lack forecast reliability assessment. A vertical field dynamic projection module enables joint modeling of ground-based and vertical observation data, and physical constraints are introduced to enhance optimization, improving the physical consistency of the forecast results. The output wide-area mean forecast field and covariance matrix provide prior information on uncertainties for subsequent probabilistic downscaling at local stations, supporting the expansion of wide-area forecasts from deterministic to probabilistic approaches and enhancing the physical reliability of the forecast results.

[0015] It should be understood that the description in this section is not intended to identify key or essential features of the embodiments of the present invention, nor is it intended to limit the scope of the invention. Other features of the invention will become readily apparent from the following description. Attached Figure Description

[0016] To more clearly illustrate the technical solutions in the embodiments of the present invention, the accompanying drawings used in the description of the embodiments will be briefly introduced below. Obviously, the accompanying drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0017] Figure 1 This is a flowchart of a multimodal heterogeneous data fusion method provided in Embodiment 1 of the present invention;

[0018] Figure 2 This is a flowchart of a wide-area atmospheric environmental factor prediction method based on multimodal fusion provided in Embodiment 2 of the present invention;

[0019] Figure 3 A flowchart of a dynamic spatial context encoding method provided in Embodiment 3 of the present invention;

[0020] Figure 4 This is a flowchart of a local site hybrid expert time series prediction method provided in Embodiment 4 of the present invention;

[0021] Figure 5 A flowchart of a probabilistic forecasting method for wide-area atmospheric pollutant concentrations provided in Embodiment 5 of the present invention;

[0022] Figure 6 This is a flowchart of a site-level atmospheric environmental factor concentration probability prediction method provided in Embodiment Six of the present invention. Detailed Implementation

[0023] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0024] Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this invention pertains. The terminology used herein in the description of this invention is for the purpose of describing particular embodiments only and is not intended to be limiting of the invention. The term "and / or" as used herein includes any and all combinations of one or more of the associated listed items.

[0025] It should be noted that some exemplary embodiments are described as processes or methods depicted as flowcharts. Although the flowcharts describe the steps as sequential processes, many of these steps can be performed in parallel, concurrently, or simultaneously. Furthermore, the order of the steps can be rearranged. A process can be terminated when its operation is complete, but it may also have additional steps not included in the figures. A process can correspond to a method, function, procedure, subroutine, subroutine, etc.

[0026] This invention provides a method for predicting atmospheric environmental factors, implemented collaboratively by a wide-area model and a local model. The wide-area model generates a broad-area coarse-resolution background field covering the target area, characterizing the regional transport and evolution of pollutants. The local model, based on the wide-area background field, integrates refined observation data from monitoring stations to output high-precision station-level predictions and probabilistic forecasts for arbitrary locations. Together, they constitute a dual-expert fusion prediction architecture of "wide-area background transport + local fine-tuning".

[0027] Based on the above architecture, the present invention illustrates the specific implementation of each technical solution through the following six embodiments. These embodiments are both independent and interconnected, collectively constituting the complete technical system of the present invention:

[0028] Example 1 relates to a multimodal heterogeneous data fusion method, which is used to fuse point spatial data and field spatial data end-to-end to generate a unified cross-modal spatiotemporal feature representation, providing high-quality input for wide-area models.

[0029] Example 2 relates to a multimodal fusion method for predicting wide-area atmospheric environmental factors, which is used to realize the background field generation function of a wide-area model and output the atmospheric environmental factor concentration forecast field at future time through a three-dimensional spatiotemporal attention encoder.

[0030] Example 3 relates to a dynamic spatial context encoding method for local site prediction, which captures the spatial correlation between sites that changes in real time with wind direction and wind speed, providing high-quality spatial features for time series prediction.

[0031] Example 4 relates to a hybrid expert time series prediction method applied to local site prediction. It uses a pluggable hybrid expert architecture to perform specialized modeling of pollution patterns at different sites and outputs site-level refined prediction results.

[0032] Example 5 relates to a probabilistic forecasting and covariance modeling method for wide-area atmospheric environmental factors, which is used to output the forecast mean field and covariance matrix in a wide-area model to achieve uncertainty quantification.

[0033] Example 6 relates to a site-level atmospheric environmental factor concentration probability prediction method, which is used to output the probability distribution parameters of any prediction location in the local model and construct the complete probability density function of atmospheric environmental factors.

[0034] For ease of understanding, the following provides a unified explanation of some of the terms used in this invention:

[0035] Spatial point data: refers to atmospheric environment-related data that exists in the form of discrete points, including but not limited to observation data from environmental monitoring stations, observation data from meteorological monitoring stations, and static information of points of interest.

[0036] Field spatial data: refers to atmospheric environment-related data that exists in the form of regular grids, including but not limited to meteorological reanalysis fields, atmospheric composition analysis fields, remote sensing observation data, and static topographic fields.

[0037] Wide-area model: refers to a model used to generate a wide-area coarse-resolution background field covering the target area, including the deterministic forecast model described in Example 2 and the probabilistic forecast model described in Example 5.

[0038] Local model: refers to a model used to output refined predictions by integrating site observation data on the basis of a wide-area background field, including the models described in Examples 3, 4 and 6.

[0039] Three-dimensional spatiotemporal attention mechanism: refers to the mechanism that captures the dynamic correlation of the atmosphere in three dimensions—time, spatial height, and spatial width—through three-dimensional window partitioning and self-attention computation.

[0040] Covariance matrix: A matrix that describes forecast uncertainty. Its diagonal elements represent the prediction variance of each factor at each location, and the off-diagonal elements represent the correlation of prediction errors between different locations or different factors.

[0041] Parametric networks: These are neural networks used to map input features to probability distribution parameters.

[0042] End-to-end training: Treat all modules of the model as a whole, drive all learnable parameters to be updated synchronously through a unified loss function, and achieve joint optimization by using gradient backpropagation.

[0043] Attention mechanism: By calculating the similarity between the query vector and the key vector, weights are generated and the value vectors are aggregated in a weighted manner to achieve adaptive filtering and fusion of information, including self-attention and cross-attention.

[0044] Residual connections: Construct parallel structures for identity mapping branches and nonlinear transformation branches to alleviate the gradient vanishing problem in deep networks. They can be combined with layer normalization to improve training stability.

[0045] Layer normalization: The mean and variance of each sample's feature dimensions are calculated independently and standardized to adapt to the feature distribution characteristics of each sample.

[0046] Uncertainty estimation refers to the process and output of characterizing the deviation between predicted atmospheric environmental factor concentrations and actual values ​​by quantifying the reliability, fluctuation range, and spatial / factor correlation of prediction results. In this invention, uncertainty estimation is achieved through a set of probability distribution parameters (such as mean and standard deviation) or a covariance matrix. Specifically, for site-level predictions, the output shows the parameters of the probability distribution (such as the mean and standard deviation of a normal distribution) followed by the atmospheric environmental factor concentration at a single prediction location, reflecting the fluctuation range and reliability of the predicted value at that location. For wide-scale predictions, the output shows a covariance matrix, which includes both the variance of the predicted values ​​at each spatial location (the magnitude of local uncertainty) and the correlation of prediction errors between different locations / different factors (global uncertainty correlation). Its core purpose is to provide a basis for reliability assessment of prediction results and support risk-sensitive environmental decision-making.

[0047] It should be noted that the "probability" in this invention refers to an objective statistical quantity learned from historical observation data based on the maximum likelihood estimation framework. This quantity is used to quantify the inherent randomness of atmospheric environmental factors and the uncertainty of model predictions. This invention directly outputs probability distribution parameters through a parameterized network and performs end-to-end training using negative log-likelihood as the loss function. This allows the model to simultaneously learn the mean and standard deviation parameters of the predictions, thereby achieving ex-post quantification of the uncertainty of the prediction results.

[0048] The above definitions are consistent throughout this invention and will not be repeated in subsequent embodiments. The specific implementations of each embodiment are described in detail below with reference to the accompanying drawings.

[0049] Example 1

[0050] This embodiment relates to a multimodal heterogeneous data fusion method, belonging to the multimodal encoding link of the wide-area model in the dual-expert fusion prediction architecture of this invention. It is used to perform end-to-end fusion of point spatial data and field spatial data to generate a unified cross-modal spatiotemporal fusion feature expression. Embodiment 1, as the encoding link in the wide-area atmospheric environmental factor prediction model (Embodiment 2), is the first step in the latter's neural network structure, used to extract, compress, and combine multimodal signals. In wide-area atmospheric environmental factor prediction, the input data usually comes from multiple heterogeneous sources, mainly including two categories: one is point spatial data, such as observation data from environmental monitoring stations, meteorological monitoring stations, and static information of points of interest; the other is field spatial data, such as remote sensing low-frequency data fields, atmospheric three-dimensional environmental observation data, and topographic static fields. These data differ significantly in spatiotemporal resolution, data form, physical meaning, and update frequency, and often exhibit time asynchrony. Existing technologies often use simple interpolation or splicing methods for processing, which makes it difficult to fully exploit the complementary information between multi-source data and easily introduces preprocessing bias.

[0051] To address the aforementioned issues, this embodiment provides a multimodal heterogeneous data fusion method. Its core concept lies in introducing an asynchronous spatiotemporal representation learning mechanism. Differentiated encoding paths are designed for the inherent characteristics of point space data and field space data. Learnable Fourier Spatiotemporal Position Encoding (FPE) is used for point-by-point / block-by-block encoding respectively. Feature compression is performed through a fixed-length compression encoder based on an attention mechanism. Information exchange is then achieved through a cross-modal fusion layer based on direct concatenation, thereby generating a fused feature representation while preserving the original structural information of each modality. The extraction method for the cross-modal spatiotemporal fusion feature representation follows a neural network model.

[0052] The specific implementation steps of this embodiment will be described in detail below with reference to the accompanying drawings.

[0053] This embodiment provides a multimodal heterogeneous data fusion method, which is implemented through a multimodal heterogeneous data fusion model to perform end-to-end fusion of different types of point spatial data and field spatial data, generating a unified cross-modal spatiotemporal feature representation, and providing high-quality input for subsequent wide-area atmospheric environmental factor prediction.

[0054] The multimodal heterogeneous data fusion model includes a data acquisition module, a point spatial data encoder, and a field spatial data encoder. The point spatial data encoder processes sparsely distributed station-type data and specifically includes an asynchronous spatiotemporal coding module (Fourier spatiotemporal location coding), a point spatial data space partitioning module, a modal-by-modal block coding module, and a fixed-length compression encoder within each block. The field spatial data encoder processes regularly gridded field data and specifically includes a field spatial data space partitioning module, a block-by-block Fourier spatiotemporal coding module, a modal-by-modal fusion compression encoder, and a field spatial data spatiotemporal location coding module. Furthermore, the model also includes a cross-modal fusion layer (stitching + MLP).

[0055] like Figure 1 As shown, the method includes the following steps:

[0056] S100, acquire multimodal input data of the target area, wherein the multimodal input data includes at least point spatial data and field spatial data.

[0057] The target area refers to the geographical region that needs to be predicted, such as one or more cities, one or more provinces, or one or more key control areas (such as industrial parks, airport perimeters), or a custom rectangular grid area defined by its latitude and longitude or projected coordinate boundaries, such as the southwest corner coordinates (x... min ,y min ) and the coordinates of the northeast corner (x max ,y max ). x min and x maxThese represent the minimum and maximum longitudes of the target area, y min and y max These are the minimum and maximum latitudes of the target area, respectively.

[0058] The point spatial data refers to data existing in the form of discrete points, characterized by sparse distribution and irregular sampling, and mainly includes at least one of the following types:

[0059] Environmental monitoring station observation data: such as time series data of atmospheric environmental factors such as PM2.5, PM10, O3, NO2, SO2, and CO collected by air quality monitoring stations;

[0060] Meteorological monitoring station observation data: such as time-series data of meteorological elements such as temperature, humidity, air pressure, wind speed, and wind direction collected by meteorological stations;

[0061] Points of interest static information: such as the geographical location of the monitoring station, the surrounding land use type, the emission source inventory, and other attribute information that does not change rapidly over time.

[0062] The field space data refers to data existing in the form of a regular grid, characterized by continuous coverage and high-dimensional density, and mainly includes at least one of the following types:

[0063] Remote sensing low-frequency data field data: such as aerosol optical thickness and vegetation index retrieved from satellites;

[0064] Three-dimensional atmospheric environmental observation data: such as three-dimensional meteorological fields and atmospheric composition fields provided by reanalysis data, including temperature, humidity, wind speed and concentration of atmospheric environmental factors at different pressure layers;

[0065] Static topographic field data: such as spatial distribution data of elevation, slope, land use type, etc.

[0066] The aforementioned multimodal input data is acquired from external data sources such as environmental monitoring centers, meteorological departments, and remote sensing data distribution platforms through the data acquisition module. After receiving the raw data, the data acquisition module also performs routine preprocessing on data with missing or abnormal values, including missing value imputation, outlier removal, and normalization, to ensure the quality and consistency of the input data and provide standardized input for subsequent encoding processing.

[0067] S110, the point space data is encoded to generate a latent space feature representation of the point space.

[0068] This step utilizes a point spatial data encoder to perform asynchronous spatiotemporal representation learning on different types of data in the point spatial data, specifically including the following sub-steps:

[0069] S1101, Fourier spatiotemporal location encoding is performed on each observation point in the point spatial data to obtain a point-by-point spatiotemporal encoding vector. The point spatial data consists of several observation points, each corresponding to a single sampling record at a certain moment, including spatial coordinates, timestamp, and observation value. Based on temporal characteristics, these observation points can be divided into two categories:

[0070] Time-series dynamic data: Observational data that changes over time, including observational data from environmental monitoring stations (such as time-series values ​​of atmospheric environmental factors such as PM2.5 and PM10 concentrations) and observational data from meteorological monitoring stations (such as time-series values ​​of meteorological elements such as temperature, humidity, and wind speed). This type of data inherently has a defined time step.

[0071] Static information data: Attribute information that does not change rapidly over time, including static information about points of interest (such as site geographic location, land use type, emission source inventory, etc.). In this embodiment, static information data does not require a separate time-series alignment module, but is assigned a unified timestamp default value (such as t=0), and then uses the same encoding method as dynamic data.

[0072] This embodiment uses learnable Fourier Spatiotemporal Position Encoding (FPE) instead of traditional position encoding. The basic idea is to first perform Fourier feature mapping on the spatial coordinates and timestamp of the observation point, then concatenate the mapped features, and finally project them onto a unified latent space through a multilayer perceptron (MLP). Specifically, for any observation point (denoted as the i-th observation point), its spatial coordinates are (x... i ,y i ,z i ), timestamp t i (This can be at any asynchronous time), the encoding vector Epos(i) is generated as follows:

[0073] Spatial encoding: Transforming the spatial coordinate vector p=[x i ,y i ,z i ] T Through a learnable frequency matrix W φ Perform a linear transformation, then calculate the sine and cosine values ​​respectively to obtain the spatial encoding vector φ(x,y,z)=Concat[sin(W φ p),cos(W φ p).

[0074] Time encoding: Encoding the timestamp t i Through another learnable frequency matrix W τ The transformation is performed and the sine and cosine values ​​are taken, while two periodic embedding features are superimposed: the day of the year doy(t) i ) and intraday hourly hod(t i), thus obtaining the time-coded vector τ(t) = Concat[sin(W τ t i ),cos(W τ t i ),Embed(doy(t i Embed(hod(t)) i ))).

[0075] Fusion and Projection: The spatial code φ and the temporal code τ are concatenated and input into a multilayer perceptron (MLP). The final output is a Fourier spatiotemporal code vector of dimension D.

[0076] Epos(i)=MLP(Concat[φ(x i ,y i ,z i ),τ(t i )])∈R D .

[0077] The parameters are defined as follows:

[0078] D: Unified latent space dimension, which is a model hyperparameter.

[0079] W φ ∈R (D / 2)×3 The learnable spatial frequency matrix maps 3D spatial coordinates to D / 2D, and then obtains D-dimensional spatial encoding through sine and cosine transformation.

[0080] W τ : A learnable time-frequency matrix, the dimensions of which are related to the time coding design.

[0081] doy(t): Converts the timestamp t to the day of the year (1~365 / 366), used to capture seasonal cycles.

[0082] hod(t): Converts the timestamp t to the hour of the day (0~23), used to capture the day-night cycle.

[0083] Embed(·): A learnable embedding layer that maps discrete date / hour values ​​to continuous vectors.

[0084] For static information data (which does not inherently have a time dimension), the same encoding method is used, but t is uniformly set. i =0 (or other fixed values). The distinction between static and dynamic data is achieved through subsequent modal embedding vectors, rather than relying on the timestamp itself.

[0085] After all the spatial observation points are encoded using the Fourier spatiotemporal position coding described above, a set of point-by-point encoded vectors carrying asynchronous spatiotemporal information is obtained.

[0086] S1102, the point-by-point encoded vector is spatially divided to obtain point space data blocks.

[0087] The target area is uniformly divided into several regular grid cells. Each observation point is assigned to a corresponding grid cell according to its spatial coordinates, and all observation points within the same grid cell constitute a point spatial data block (Patch).

[0088] In the spatial partitioning (grid mapping) operation of this step, the target area is uniformly divided into Gw columns along the longitude (x) direction and Gh rows along the latitude (y) direction, resulting in a total of Gh×Gw regular grid cells. Each grid cell corresponds to a point spatial data patch. For each coded observation point, based on its coordinates (x... i y i Calculate the row index r of the corresponding grid cell. i and column index c i :

[0089] r i =min(roundup((y i -y min ) / △y), Gh-1), c i =min(roundup((x i -x min ) / △x), Gw-1), roundup() means rounding up, △x and △y represent the spatial step size of the grid cell, Δx=(x max -x min ) / Gw, Δy=(y max -y min ) / Gh.

[0090] The global patch index of this observation point is p. i =r i ×Gw+c i All observations assigned to the same patch constitute the input matrix of that patch.

[0091] The mesh size Gh, Gw must be aligned with the total number of blocks in the field space data. See S1201 for specific alignment methods.

[0092] If there are no observation points in a certain grid cell (empty patch), bilinear interpolation is used to fill in the gaps from the observation points of the non-empty patches in the 8-neighborhood; if there are no data in the neighborhood, the global mean of the mode is used to fill in the gaps.

[0093] S1103, perform modal projection and modal embedding on each point space data block to obtain the corresponding point space input sequence.

[0094] For the j-th point space data block (j ranges from 1 to Np1, where Np1 is the total number of point space data blocks), it contains N j There are 10 observation points. Each observation point carries a Fourier spatiotemporal code Epos(i) and the original feature vector.

[0095] First, the original features of each observation point are mapped to the latent space D through a modality-specific linear projection layer, resulting in a corresponding projection vector. Then, this projection vector is added to the Fourier spatiotemporal code Epos(i) to obtain a fusion vector. Finally, to distinguish between different modalities (e.g., ground station and mobile observation vehicle), a learnable modality identifier vector is superimposed on the fusion vector to obtain the final vector. The final vectors of all observation points within the j-th spatial data block are arranged in spatial order to obtain the input sequence Z. in (j) serves as the input to the subsequent compression encoder.

[0096] S1104, perform feature compression on each point space data block to obtain the point space latent space feature representation.

[0097] The feature compression includes: dividing the encoded feature sequence into data blocks according to spatial location, and compressing the variable-length sequence in each data block into a predetermined number of latent vectors using an attention-based compression encoder. The compression encoder uses a set of learnable latent query vectors as queries, and the feature sequences within the data block as keys and values. It compresses the variable-length input sequence into intermediate latent vectors using a cross-attention mechanism, and applies a self-attention mechanism to the intermediate latent vectors for feature enhancement to obtain the predetermined number of latent vectors.

[0098] Specifically, the input sequence Z of each data block in (j) Input block fixed-length compression encoder. This encoder employs a perceptual resampler architecture based on a cross-attention mechanism: with a set of learnable latent query variables Z lat ∈R k×D (k is a fixed number, and k is much smaller than N) j ) as the query, with Z in (j) serves as the key and value. Cross-attention is used to compress the variable-length sequence into a fixed-length output, yielding latent vectors. Self-attention is then applied to these latent vectors for feature enhancement. Ultimately, each data block outputs a dimension of k×D, meaning each data block is compressed into k fixed-dimensional latent vectors. The outputs of all data blocks are stacked to obtain the point-space latent feature representation Spoint∈R. Np1×k×D The encoder supports asynchronous incremental updates: when new observations arrive, forward propagation continues with the current latent variables as the initial state.

[0099] Through the above processing, each point spatial data block is compressed into k latent vectors, and the latent vectors of all spatial blocks together constitute the latent space feature representation of the point space. This feature representation not only preserves the local spatial structure of the original point data, but also removes redundant information through compression, providing efficient feature input for subsequent cross-modal fusion.

[0100] S120, the field space data is encoded to generate a latent space feature representation of the field space.

[0101] This step utilizes a field space data encoder to extract features from the regularly gridded field data, specifically including the following sub-steps:

[0102] S1201, The field space data is spatially divided to obtain field space data blocks.

[0103] The field space data (size H×W, sharing a single timestamp t) covering the target area is directly divided into non-overlapping S×S pixel blocks along the spatial dimension, with a total number of blocks Np = H / S×W / S. To align with the number of patches in the point space data blocks, the spatial grid size for point space should be set to Gh = H / S and Gw = W / S, ensuring that the number of blocks in both parts is equal. For 3D field data containing a vertical dimension, blocks are cut along the (H,W) plane, preserving the vertical dimension as a channel and flattening it. After field data of different resolutions are independently cut into blocks, the number of patches is aligned by interpolation or adaptive step size adjustment.

[0104] S1202 performs block-by-block Fourier spatiotemporal coding on each field space data block.

[0105] For the p-th data block (p=1,2,…,Np2, where Np2 is the total number of field space data blocks), calculate its spatial barycenter coordinates (x...). p ,y p (For example, take the average of the coordinates of all grid points within the block), and using the shared timestamp t of the field data, generate the encoding vector Epos(p)∈R using the same Fourier spatiotemporal encoding function as the point space. D It should be noted that the Fourier encoding of the field space is performed after the blocks are sliced, because the entire field has only one timestamp and does not need to be encoded point by point.

[0106] S1203, Modal projection is performed on each field space data block to obtain the corresponding field space input sequence.

[0107] For the p-th data block, it contains S 2 Each pixel / grid point carries a Fourier spatiotemporal code Epos(p) (obtained from step S1202) and the original feature vector.

[0108] First, the original features of each pixel / grid point are mapped to the latent space D through a mode-specific linear projection layer, obtaining the corresponding projection vector. Then, this projection vector is added to the Fourier spatiotemporal code Epos(p) to obtain the fused vector. Finally, to distinguish different field spatial modes (such as satellite cloud images and reanalysis fields), a learnable mode identifier vector is superimposed on the fused vector to obtain the final vector. The final vectors of all pixels / grid points within the p-th field spatial data block are arranged in spatial order to obtain the input sequence Z. in (p), which serves as the input for the subsequent compression encoder.

[0109] S1204, perform feature compression on each field space data block to obtain the latent space feature representation of the field space.

[0110] The input sequence of each field space data block is fed into a modal fusion compression encoder. This encoder has the same architecture as the point space compression encoder (but with independent parameters), and also employs a cross-attention-based perceptual resampler architecture to learnable latent variables from a fixed-length S... 2 The input is compressed into a fixed-dimensional k×D output. The outputs of all data blocks are stacked to obtain the field space latent space feature representation Sfield∈R. Np×k×D .

[0111] S130, the latent space feature representation of the point space and the latent space feature representation of the field space are fused to generate a cross-modal spatiotemporal fusion feature representation.

[0112] This step fuses the latent space features output by the point space encoder and the field space encoder through a cross-modal fusion layer. Since S1102 and S1201 have ensured that the total number of data blocks in the point space and the field space are the same, and that each data block index corresponds to the same spatial location, features from both in the same spatial region can be fused directly.

[0113] The fusion strategy employs direct stitching combined with a multilayer perceptron (MLP), rather than a cross-attention mechanism. The specific processing flow is as follows:

[0114] For each data block, extract the latent space features in the point space and the latent space features in the field space corresponding to that block. Both have the same dimensional structure: each data block contains a fixed number of k latent vectors, and each latent vector has a dimension of D. These two sets of features are concatenated along the feature dimensions of the latent vectors, that is, each latent vector is expanded from the original D dimensions to 2D dimensions, but the number of latent vectors k remains unchanged.

[0115] Next, the concatenated features are input into a multilayer perceptron (MLP). This MLP contains at least one hidden layer with an output dimension of D, which projects the concatenated features back into the original unified latent space dimension D. The MLP can also include conventional components such as batch normalization, activation functions (such as GELU), and Dropout to enhance expressive power and training stability.

[0116] After performing the above operations on all data blocks sequentially, they are stacked in spatial order to form a unified spatiotemporal bottleneck tensor. The size of this tensor is: total number of data blocks × number of latent vectors k in each data block × unified feature dimension D. In this tensor, the k latent vectors corresponding to each data block encode the point-source-field-source joint features of that spatial location, and these latent vectors have established internal correlations through the self-attention layer of the preceding perceptual resampler. This tensor can be directly connected to downstream prediction heads (e.g., decoders based on Transformer or convolutional structures) for meteorological element forecasting.

[0117] Compared to fusion methods based on cross-attention, this solution uses direct splicing plus MLP, which has the following advantages:

[0118] Lower computational complexity: The computational cost of cross-attention increases with the square of the product of the number of data blocks and the number of latent vectors, while the computational cost of direct concatenation increases only linearly with the product of the number of data blocks, the number of latent vectors, and the feature dimension, significantly reducing computational overhead.

[0119] Avoid modal competition: Cross attention may cause point space features and field space features to inhibit each other (i.e., one modality dominates and the information of the other modality is ignored), while direct splicing preserves the original information of both, and the MLP learns how to combine them optimally, which is more stable and reliable.

[0120] Simple to implement and easy to train: No need to maintain cross-modal attention weight matrix, reducing the number of parameters and the risk of overfitting.

[0121] Through the above fusion process, the local refined observation information in the point spatial data is injected into the field spatial data, realizing the complementary advantages of the two modes and providing high-quality input features for subsequent prediction tasks.

[0122] In this embodiment, the temporal information of both point space data and field space data has been embedded separately through the preceding Fourier spatiotemporal location encoding steps (S1101 and S1202). Specifically, the independent timestamp of each observation point in the point space has been injected during point-by-point encoding; the shared timestamp of the entire field in the field space has been injected during block-by-block encoding. Therefore, this method does not require adding an additional time encoding layer before fusion, nor does it require forced time alignment of data from different sources. The above design enables the multimodal heterogeneous data fusion model proposed in this invention to be natively compatible with asynchronous observation data commonly found in actual business operations.

[0123] Furthermore, in this embodiment, feature compression supports asynchronous incremental updates: when new observation data is acquired, the query variables of the current compression encoder remain unchanged, and forward propagation is performed only on the new data to update the latent space feature representation of the corresponding path. In this way, new observation information can be efficiently integrated into existing fused features without reprocessing all historical data, reducing redundant calculations and making it suitable for real-time or streaming data scenarios.

[0124] Furthermore, the method employs an end-to-end training approach, where learnable parameters in all steps are jointly optimized through gradient backpropagation.

[0125] End-to-end training refers to treating the entire multimodal heterogeneous data fusion model (including the point spatial encoder, field spatial encoder, cross-modal fusion layer, all embedding layers, attention parameters, etc.) as a whole, and using a unified loss function to drive the synchronous update of all learnable parameters. Specifically:

[0126] The range of learnable parameters includes, but is not limited to, the spatial frequency matrix and temporal frequency matrix in Fourier spatiotemporal location coding, the linear projection layer weights specific to each modality, the cross-attention projection matrix in the perceptual resampler of point space and field space, as well as learnable latent variables, self-attention layer parameters, modality identifier vectors, multilayer perceptron (MLP) weights in cross-modal fusion layers, and other optional layer normalization parameters.

[0127] Joint optimization process:

[0128] Forward propagation: Multimodal asynchronous data from the training samples are sequentially input into the model, passing through a point spatial encoder, a field spatial encoder, and a cross-modal fusion layer, ultimately outputting a cross-modal spatiotemporal fusion feature representation. This feature representation is then input into downstream prediction task modules (such as meteorological element prediction heads) to obtain prediction results (such as the PM2.5 concentration distribution at a future time).

[0129] Loss Calculation: A loss function is constructed based on the difference between the predicted results and the actual observed values. For example, mean squared error (MSE) can be used for regression tasks, while cross-entropy loss can be used for classification tasks. Multiple loss terms can be stacked as needed (such as spatial weighting, multi-time-step output, etc.).

[0130] Backpropagation: Calculates the gradient of the loss function with respect to each learnable parameter in the model using automatic differentiation techniques.

[0131] Parameter update: Use gradient descent optimizers (such as Adam and SGD) to update the parameters of all modules at once.

[0132] Technical benefits: End-to-end joint optimization avoids the error accumulation and local optima problems that may occur in traditional staged training. Driven by a unified goal, each module adjusts collaboratively, learning the optimal intra-modal feature extraction method and inter-modal interaction weights, thereby generating higher-quality fused feature representations and improving the accuracy of the final prediction task.

[0133] Through the end-to-end training method described above, all learnable parameters in this method can be jointly optimized, ensuring the optimal overall performance of the model.

[0134] The multimodal heterogeneous data fusion method provided in this embodiment designs differentiated encoding paths for point space data and field space data, introduces asynchronous spatiotemporal representation learning and Fourier spatiotemporal position encoding, and employs a fixed-length compression encoder based on cross-attention for feature compression. Finally, it achieves low-complexity merging through direct concatenation and fusion. This method has the following technical features:

[0135] (1) Preserve the original observation fidelity and reduce information loss. Point spatial data adopts Fourier encoding with independent timestamps for each point, which does not require forced alignment to a unified time axis through interpolation, thus avoiding the smoothing of extreme weather event signals to a certain extent. Field spatial data retains the original grid resolution and is only locally organized by slicing.

[0136] (2) Supports asynchronous time data. Fourier spatiotemporal location coding maps irregular acquisition times to a high-dimensional orthogonal basis space. The model can calculate the spatiotemporal similarity between any two asynchronous observations, achieving "fuzzy matching" without the need for forced time alignment of the input data.

[0137] (3) Reduce computational complexity. The fixed-length compression encoder compresses the input sequence into a fixed dimension through learnable latent variables, reducing the computational complexity from the quadratic level of self-attention to a linear relationship with the input length; the fusion layer uses direct splicing instead of cross-attention, avoiding secondary complexity.

[0138] (4) Supports flexible combination of multimodal data. Through a unified block size and modal projection layer, multimodal data from different sources and with different resolutions can be combined. The compression encoders in point space and field space are executed independently. When adding a new modality, only the corresponding encoder and modal embedding need to be added.

[0139] (5) Supports streaming asynchronous updates. During inference, when new observation data arrives (such as delayed reporting from a ground station), only one forward propagation needs to be performed on the compressed encoder of that mode to update the corresponding latent variables. The fusion layer can then recalculate to generate new predictions without recalculating the entire model.

[0140] (6) It possesses the ability to tolerate missing data. The point space can be filled in by neighborhood interpolation or global mean, and the continuous coverage characteristic of the field space can be used to supplement the sparse point space across modes. The model can maintain performance even when some sensors are missing or delayed. In summary, this embodiment provides input feature representations for atmospheric environmental factor prediction tasks through asynchronous spatiotemporal representation learning, Fourier spatiotemporal location coding, fixed-length compressed encoder, and low-complexity direct splicing fusion.

[0141] Based on the same inventive concept, embodiments of the present invention also provide a multimodal heterogeneous data fusion apparatus, comprising:

[0142] The input module is used to acquire multimodal input data of the target area, wherein the multimodal input data includes at least point spatial data and field spatial data;

[0143] The point space encoding module is used to perform asynchronous spatiotemporal encoding and feature compression on the point space data to generate a latent space feature representation of the point space.

[0144] The field space coding module is used to perform synchronous spatiotemporal coding and feature compression on the field space data to generate a latent space feature representation of the field space.

[0145] The fusion module is used to fuse the latent space features of point space and field space at the same spatial location to generate a cross-modal spatiotemporal fusion feature representation.

[0146] Optionally, the device further includes at least one of the following features:

[0147] The point spatial encoding module is specifically used for: performing learnable Fourier spatiotemporal position encoding on each observation point in the point spatial data using its independent spatial coordinates and acquisition timestamp to obtain a point-by-point spatiotemporal encoding vector; mapping the point-by-point encoding vector to regular grid cells according to spatial position, with each grid cell constituting a point spatial data block; and performing fixed-length compression on the variable-length sequence within each data block using a perceptual resampler based on cross-attention to obtain a point spatial latent space feature representation.

[0148] The field space encoding module is specifically used to: divide the field space data into regular data blocks along the spatial dimension; perform learnable Fourier spatiotemporal position encoding on each data block using its spatial centroid coordinates and the timestamp shared by the field data to obtain a block-by-block spatiotemporal encoding vector; and compress the encoding vector sequence within each data block to a fixed length using a perceptual resampler based on cross-attention to obtain a latent space feature representation of the field space.

[0149] The Fourier spatiotemporal location encoding includes a learnable spatial frequency matrix and a temporal frequency matrix. The temporal encoding also includes seasonal embedding and day-night embedding. The encoding result is projected onto a unified latent space via a multilayer perceptron.

[0150] The spatial division granularity of the point spatial data block is determined by the segmentation granularity of the field spatial data block, so that the number of point spatial data blocks and field spatial data blocks are equal and their spatial positions correspond one-to-one.

[0151] The fusion module is specifically used for: for point spatial data blocks and field spatial data blocks at the same spatial location, directly concatenating their latent vector sets on the feature dimension, and then obtaining the fusion feature at that location through nonlinear mapping by a multilayer perceptron; stacking the fusion features of all spatial locations in spatial order to form a spatiotemporal bottleneck tensor.

[0152] The device employs an end-to-end training method, where the learnable parameters of all modules are jointly optimized through gradient backpropagation.

[0153] Example 2

[0154] This embodiment relates to a multimodal fusion-based wide-area atmospheric environmental factor prediction method, belonging to the wide-area modeling link in the dual-expert fusion prediction architecture of this invention. It is implemented through a multimodal fusion-based wide-area atmospheric environmental factor prediction model, used to generate future atmospheric environmental factor forecast fields at a wide-area scale, providing macroscopic background field support for regional pollution transport analysis, emergency decision-making, and local refined prediction (Examples 3 to 6). The core challenge in wide-area atmospheric environmental factor prediction lies in how to efficiently and accurately capture the complex dynamic evolution of the atmospheric system in horizontal space, vertical stratification, and time dimensions. While traditional numerical simulation methods (such as WRF-Chem and CMAQ) can generate continuous forecast fields, their computational cost is enormous, often requiring several hours to complete a single forecast, making it difficult to meet real-time requirements. Furthermore, these methods rely on manually designed physical parameterization schemes, making it difficult to adaptively learn and optimize from massive historical observation data. Existing purely data-driven methods (such as models based on convolutional neural networks) are limited by local receptive fields, making it difficult to effectively capture long-distance transport of pollutants over thousands of kilometers. Furthermore, they do not make sufficient use of vertically layered information, often simply stacking multiple layers of data into channels for processing, ignoring the physical correlation in the vertical direction, which limits the accuracy of forecasts.

[0155] The multimodal fusion-based wide-area atmospheric environmental factor prediction model comprises, in sequence, an input layer, the multimodal heterogeneous data fusion model described in Example 1, a three-dimensional spatiotemporal sensing encoder, and a decoder. The multimodal heterogeneous data fusion model described in Example 1 is used to perform end-to-end fusion of point spatial data and field spatial data to generate a cross-modal spatiotemporal fusion feature representation. The three-dimensional spatiotemporal sensing encoder is used to extract spatiotemporal features from the input data, capturing the dynamic correlation of the atmosphere in horizontal space, vertical stratification, and time dimension. The decoder is used to map the spatiotemporal feature map output by the encoder into a future atmospheric environmental factor forecast field. This model jointly optimizes the parameters of all modules through end-to-end training, achieving a direct mapping from multi-source input to wide-area forecast.

[0156] To address the aforementioned issues, this embodiment provides a multimodal fusion method for predicting wide-area atmospheric environmental factors. Its core concept lies in using regularly gridded field spatial data (or the fusion features of point spatial data and field spatial data processed in Embodiment 1) as input. A three-dimensional spatiotemporal attention encoder simultaneously captures the dynamic correlations of the atmosphere in horizontal space, vertical stratification, and time, generating a multi-scale spatiotemporal feature map. A decoder then outputs the forecast field of atmospheric environmental factors for future times. This method combines the efficiency of deep learning with the accuracy of three-dimensional spatiotemporal modeling, enabling wide-area forecasts to be completed within minutes, providing a high-quality macroscopic background field for subsequent local refined predictions.

[0157] The specific implementation steps of this embodiment will be described in detail below with reference to the accompanying drawings.

[0158] like Figure 2 As shown, the method includes the following steps:

[0159] S200, acquire multi-source input data for historical time steps, wherein the multi-source input data includes at least field space data, or includes fusion features of point space data and field space data after fusion processing.

[0160] In this embodiment, the multi-source input data includes field space data or fusion features of multiple consecutive historical time steps (e.g., the current time T and its previous time T-1), which are used to capture the short-term temporal evolution of atmospheric state.

[0161] The type and acquisition method of the multi-source input data are the same as those of the multimodal input data described in Embodiment 1. Specifically, the multi-source input data includes at least one of the following two forms:

[0162] Field spatial data: Its definition is the same as in Example 1, referring to atmospheric environment-related data that exists in the form of a regular grid. It has the characteristics of continuous coverage and high-dimensional density, including meteorological reanalysis fields, atmospheric composition analysis fields, remote sensing observation data, and static field data.

[0163] Fusion Features: These refer to the cross-modal spatiotemporal fusion feature expressions generated after processing by the multimodal heterogeneous data fusion method described in Example 1. These fusion features have undergone differentiated processing by point spatial data encoders and field spatial data encoders, enabling end-to-end fusion of sparse monitoring station observation data with dense field spatial data. This preserves both the local, accurate observation information of the station data and the spatial continuity of the field spatial data, providing more accurate input for wide-area prediction.

[0164] Multi-source input data is obtained from external data sources such as environmental monitoring centers, meteorological departments, and remote sensing data distribution platforms through the input layer of the prediction model, or directly received from the output of the fusion model in Example 1. Before being input to the 3D spatiotemporal sensing encoder, the input layer also needs to perform necessary preprocessing on the data, including spatiotemporal alignment (unifying data from different sources to the same spatial grid and temporal resolution), missing value imputation, outlier removal, and normalization, to ensure the quality and consistency of the input data and provide standardized input for subsequent spatiotemporal feature extraction steps.

[0165] It should be noted that when using only field space data as input, this method can still run independently and generate forecast results; however, when using the features fused in Example 1 as input, the forecast accuracy can be further improved due to the incorporation of real information from station observations. This flexible input design allows this method to adapt to application scenarios under different data conditions.

[0166] S210, the multi-source input data is input to the three-dimensional spatiotemporal sensing encoder, and spatiotemporal features are extracted through the three-dimensional window attention mechanism to generate a spatiotemporal feature map.

[0167] This step is the core feature extraction step for wide-area atmospheric environmental factor prediction. Its goal is to perform deep spatiotemporal coding on the multi-source input data acquired by S200, capture the complex dynamic correlation of atmospheric pollutants in horizontal space, vertical stratification and time dimension, and provide rich feature representation for subsequent forecast field generation.

[0168] (1) Overall architecture of the three-dimensional spatiotemporal sensing encoder

[0169] The 3D spatiotemporal sensing encoder is built on a 3D vision Transformer architecture and consists of multiple stacked encoding stages. Each encoding stage contains several 3D Swing Transformer blocks, and the receptive field is gradually expanded through downsampling operations to achieve multi-scale feature extraction from local to global. The encoder's input is multi-source input data output from the S200 (shape: batch size × input time step × spatial height × spatial width × number of channels), and the output is a spatiotemporal feature map with the same spatial size but encoded channel count.

[0170] Unlike conventional 2D visual Transformers, atmospheric forecast data involves complex vertical structures—the same variable may have physical relationships across different atmospheric layers (such as the boundary layer, troposphere, and stratosphere), but also carry distinct physical meanings. Traditional 2D ViT structures simply stack vertical layers as channels for processing, making it difficult to efficiently model these complex vertical relationships. Therefore, this embodiment employs an innovative 3D ViT multi-scale spatiotemporal sensing architecture, extending the spatial dimension from 2D to 3D, simultaneously processing horizontal space and vertical layers to more deeply mine information about atmospheric states at various scales.

[0171] (2) Three-dimensional window attention mechanism

[0172] To simultaneously capture atmospheric dynamics across three dimensions—temporal, spatial altitude, and spatial width—this embodiment extends the standard two-dimensional window attention to a three-dimensional window attention. Specifically:

[0173] 3D Window Partitioning: The input feature map is divided into non-overlapping 3D windows, each with dimensions Twindow × Hwindow × Wwindow, where Twindow is the temporal window size, and Hwindow and Wwindow are the window sizes in the spatial height and width directions, respectively. In this embodiment, since the number of input time steps is fixed at 2, Twindow can be set to 2, allowing each window to simultaneously contain local spatiotemporal regions of two historical time steps, facilitating the capture of short-term evolution features.

[0174] Intra-window self-attention computation: Within each 3D window, all feature points within the window are flattened into a sequence, and standard self-attention computation is performed. By introducing a relative positional bias, the model can perceive the spatiotemporal positional relationships of feature points within the window. This design retains the ability of self-attention to model long-range dependencies while reducing the computational complexity from a global O(T×H×W). 2 Reduce to window-level O((Twindow×Hwindow×Wwindow)) 2 ×T×H×W / (Twindow×Hwindow×Wwindow)), which makes it possible to process high-resolution spatiotemporal data.

[0175] Window sliding mechanism: To enhance information interaction between different windows, a window sliding operation is introduced between consecutive 3D SwinTransformer blocks. For example, the first block uses a regular window division, and the second block offsets the window in the time, height, and width directions by (Twindow / 2, Hwindow / 2, Wwindow / 2) respectively, thereby establishing a cross-window connection and realizing the spatial propagation of information.

[0176] (3) Specialized processing of vertical layering

[0177] Considering that atmospheric vertical stratification has clear physical meaning (such as boundary layer, troposphere, and stratosphere), and that the correlation between different altitude layers varies with meteorological conditions, this embodiment introduces a vertical dimension perception mechanism into the three-dimensional window attention. Specifically, one or a combination of the following two methods can be used:

[0178] Option 1: Vertical Channel Fusion. Vertical observation data is stitched together with the ground observation field along the channel dimension as an additional channel, allowing the 3D window to simultaneously process information from all vertical layers. This method is simple and direct, suitable for scenarios with a small number of vertical layers and relatively fixed inter-layer relationships.

[0179] Option 2: 3D Convolutional Preprocessing. Before inputting the data into the encoder, several 3D convolutional layers are used to extract local features from the data containing the vertical dimension, initially fusing the vertical information. This is then input into a 3D SwingTransformer block for deep modeling. This approach can explicitly capture the local structure in the vertical direction and is more suitable for scenarios with a large number of vertical layers or complex inter-layer relationships.

[0180] The preferred embodiment is the second scheme, which first maps the input tensor to a higher-dimensional feature space through a three-dimensional convolutional layer, while performing local feature aggregation in the vertical direction, and then feeds it into the subsequent three-dimensional Swing Transformer block.

[0181] (4) Multi-scale hierarchical feature extraction

[0182] To simultaneously capture multi-scale atmospheric phenomena such as local emissions, regional transport, and wide-area circulation, the encoder employs a hierarchical downsampling design. Each encoding stage consists of several 3D Swing Transformer blocks and a downsampling layer:

[0183] Phase 1: Input resolution H×W, time step Tin, number of channels C1. Fine-grained spatiotemporal features are extracted using several 3D SwingTransformer blocks.

[0184] Phase Two: After the downsampling layer, the spatial resolution is reduced to H / 2×W / 2, the time steps remain unchanged or are merged as needed, and the number of channels increases to C2. Medium-scale features are then extracted again using a 3D Swing Transformer block.

[0185] Phase 3: Further downsampling to H / 4×W / 4, increasing the number of channels to C3, to capture large-scale regional transmission characteristics.

[0186] Downsampling can be achieved using 3D convolution or 3D patch merging (e.g., linearly projecting concatenated 2×2×2 neighborhood features). Through this hierarchical structure, the spatiotemporal feature map output by the encoder encodes atmospheric state information at different scales, providing rich multi-scale features for subsequent mean field and covariance prediction.

[0187] (5) Output spatiotemporal feature map

[0188] After multiple stages of encoding, the final output is a spatiotemporal feature map, the spatial resolution of which can be selected according to task requirements (e.g., retaining the original resolution or using features from a certain intermediate layer). This spatiotemporal feature map contains deep spatiotemporal features extracted from the input spatiotemporal cube, preserving both local details and incorporating global context, and will serve as input for subsequent mean prediction heads and covariance prediction heads.

[0189] It should be noted that the hyperparameters such as the size of the three-dimensional window, the number of encoder stages, and the number of channels can be adjusted according to factors such as the actual prediction area size and computing resources. Those skilled in the art can make reasonable configurations without departing from the concept of this invention. Through the processing of the three-dimensional spatiotemporal sensing encoder, this method can effectively capture pollution transport over thousands of kilometers, vertical material exchange, and temporal evolution patterns, laying the foundation for high-precision wide-area forecasting.

[0190] S220, based on the spatiotemporal feature map, output the atmospheric environmental factor forecast field for future times through the decoder.

[0191] This step utilizes a decoder to map the spatiotemporal feature map extracted by S210 into a wide-area atmospheric environmental factor forecast field for future times, achieving decoding from feature space to physical space. The decoder adopts a multi-head joint output structure, simultaneously learning the concentration distribution of multiple atmospheric environmental factors through a shared bottom feature layer, fully utilizing the physical correlation and statistical correlation between different factors to improve overall forecast accuracy.

[0192] (1) Decoder overall architecture

[0193] The decoder is symmetrically designed with the 3D spatiotemporal-aware encoder and consists of multiple upsampling stages, progressively restoring spatial resolution and reducing the number of channels, ultimately outputting a prediction field with the same grid size as the input space. Each upsampling stage includes an upsampling layer (such as transposed convolution, bilinear interpolation, or pixel reconstruction) and several convolutional layers to fuse high-level semantic features and supplement detailed information.

[0194] (2) Multi-head joint output structure

[0195] To simultaneously predict multiple atmospheric environmental factors, the decoder employs a multi-head joint output structure in its final layer. Specifically, the decoder's shared bottom feature layer first progressively upsamples the spatiotemporal feature map to the original input resolution, resulting in a shared feature map of shape H×W×D (where D is the shared feature dimension). Then, for each atmospheric environmental factor, an independent output head (i.e., a convolutional layer) is set up to map the shared feature map to the concentration prediction field of that factor.

[0196] Yc=Conv c (F shared ), c=1,2,…,C.

[0197] Where C represents the number of predicted atmospheric environmental factors, Conv c A dedicated convolutional layer (typically a 1×1 or 3×3 convolution) is used for the c-th factor, outputting a single-channel forecast field with shape H×W×1. The forecast fields from all output heads are concatenated along the channel dimension to form the final multi-factor forecast field Y∈R. H×W×C .

[0198] The atmospheric environmental factors include at least two or more of PM2.5, PM10, O3, NO2, SO2, and CO. In practical applications, the types of factors to be predicted can be flexibly selected according to the task requirements.

[0199] (3) Advantages of multi-task learning

[0200] By sharing a low-level feature layer and simultaneously learning multiple factors, the decoder can automatically uncover implicit correlations between different pollutants (e.g., the correlation between PM2.5 and PM10, and the chemical coupling between O3 and NO2), thus forming a synergistic learning effect during training. This multi-task joint learning mechanism has the following advantages compared to building a separate model for each factor:

[0201] High parameter efficiency: The shared feature layer significantly reduces the total number of parameters, thus reducing the risk of overfitting;

[0202] High computational efficiency: All factors can be predicted in a single forward propagation, resulting in fast inference speed;

[0203] Improved accuracy: The gradient information of different factors complements each other, which helps to learn more robust feature representations, especially in factors or regions with sparse data.

[0204] (4) Output post-processing and verification

[0205] The forecast field output by the decoder represents the concentration distribution of atmospheric environmental factors at future times. To ensure the physical validity of the forecast results, necessary post-processing can be performed on the output, such as spatial smoothing and threshold pruning (to ensure non-negative concentrations). Furthermore, a validation mechanism can be introduced during training to periodically evaluate the forecast accuracy of each factor and adjust the model structure or hyperparameters as needed.

[0206] Through the above steps, this embodiment realizes end-to-end mapping from multi-source input to a wide-area multi-factor forecast field, providing a high-quality macro background field for subsequent environmental management decisions and local refined forecasting.

[0207] Furthermore, the method also includes a parameter fine-tuning step: when optimization is needed for a specific environmental factor or a specific region, the parameters of the three-dimensional spatiotemporal sensing encoder are frozen, and only the output head of the decoder corresponding to the specific factor is fine-tuned.

[0208] In practical applications of wide-area atmospheric environmental factor prediction, the following two typical optimization requirements may exist:

[0209] Optimization for specific environmental factors: For example, in a region where ozone pollution is prominent, it is necessary to improve the model's prediction accuracy for O3, while the accuracy requirements for other factors are relatively lower.

[0210] Optimization for specific regions: For example, when transferring a pre-trained model to a new target region (such as from East China to South China), it is necessary to quickly adapt to the geographical features and emission characteristics of the new region.

[0211] To meet the above requirements, this embodiment provides a parameter fine-tuning strategy. When optimization is needed for a specific environmental factor or a specific region, the following steps are performed:

[0212] Freeze the backbone network: Keep all parameters of the 3D spatiotemporal sensing encoder unchanged, that is, freeze the learned general spatiotemporal feature extraction capabilities. These parameters have already grasped the basic laws of atmospheric dynamic evolution in large-scale data pre-training, have good generalization ability, and do not need to be relearned for specific tasks.

[0213] Fine-tuning the decoder output head: Only the output head corresponding to a specific environment factor in the decoder is updated. Specifically, in a multi-head joint output structure, each environment factor corresponds to an independent output head (usually a 1×1 convolutional layer or a lightweight network layer). Through the backpropagation algorithm, only the output head parameters of the target factor are updated, while the output head parameters of other factors remain unchanged.

[0214] Lightweight training: Lightweight training is performed using a small amount of labeled data for the target region or target factor. Because the number of trainable parameters is greatly reduced (only 1% to 5% of the original model), the amount of training data and computational resources required are also reduced accordingly. Fine-tuning can usually be completed within a few hours, and it is less prone to overfitting.

[0215] Optional full fine-tuning: If the performance still does not meet the requirements after fine-tuning, some encoder layers (such as the last few layers) can be unfrozen for joint optimization. However, this embodiment recommends fine-tuning only the output head to retain the universality of the pre-trained knowledge to the greatest extent.

[0216] This efficient parameter fine-tuning strategy has the following technical advantages:

[0217] Low optimization cost: It can be completed with only a small amount of data and computing resources, making it suitable for frequent model update needs in actual business.

[0218] Avoid destructive forgetting: The parameters of the backbone network remain unchanged, so the learned general spatiotemporal features are not lost;

[0219] Fast response: It can quickly optimize specific factors for newly emerging pollution problems (such as sudden sandstorms);

[0220] Multi-scenario adaptation: The same pre-trained model can be fine-tuned by different output heads to serve multiple downstream tasks (such as forecasting services in different regions) at the same time, reducing deployment costs.

[0221] It should be noted that this strategy applies to the deterministic forecasting model described in Example 2 (only the mean prediction head corresponding to a specific factor in the decoder is fine-tuned), and also to the probabilistic forecasting model described in Example 5 (the mean prediction head or covariance prediction head can be fine-tuned, selected according to optimization needs). Those skilled in the art can select the specific modules and number of layers for fine-tuning according to actual needs, all of which fall within the protection scope of this invention.

[0222] Furthermore, the method also includes a physical principle alignment step: introducing a physical-guided machine learning paradigm during model training, and constraining the deep learning model by combining atmospheric physical laws or intermediate information from numerical simulation models, so that the prediction results follow atmospheric physical and chemical principles.

[0223] During model training, physical constraint terms can be introduced into the loss function to penalize predictions that violate physical laws in a differentiable form. These physical constraints can include at least one of the following: mass conservation constraints, chemical kinetic constraints, and boundary layer mixing constraints. For a detailed description of the implementation, please refer to the "Physical Constraint Enhancement Optimization Steps" in Example 5; here, only a conceptual explanation is provided: mass conservation constraints constrain the divergence of pollutant flux based on the continuity equation; chemical kinetic constraints apply empirical relationships to the photochemical reaction system; and boundary layer mixing constraints adjust the vertical concentration gradient. These physical constraint terms are added to the total loss function in a weighted form, guiding the model parameters to update in a direction that satisfies physical laws through backpropagation. This improves the physical reliability of the prediction results while maintaining the efficiency of data-driven approaches.

[0224] Furthermore, the method also includes an uncertainty quantification step: the covariance matrix of the forecast uncertainty is output by the covariance prediction head, and together with the mean field output by the decoder, it defines the multivariate Gaussian distribution of the atmospheric state at future times. In wide-area atmospheric environmental factor prediction, simply outputting a deterministic concentration mean field is insufficient to meet the requirements for evaluating the reliability of forecast results in practical applications. Due to the inherent errors in the meteorological input field, the bias in model approximations, and the inherent chaotic characteristics of the atmospheric system, forecast results at different spatial locations and at different times have different degrees of uncertainty, and the prediction errors between different locations often exhibit spatial correlation. Therefore, this embodiment adds a covariance prediction head to the mean field output by the decoder to quantify the forecast uncertainty.

[0225] (a) Structure of the covariance prediction head

[0226] The covariance prediction head is connected in parallel with the decoder, taking the spatiotemporal feature map extracted by S210 as input, but outputting covariance matrix parameters describing the prediction uncertainty. The covariance prediction head can adopt a convolutional network structure similar to the mean prediction head, mapping the high-dimensional spatiotemporal feature map into covariance parameters through several convolutional layers.

[0227] (II) Mathematical meaning of the covariance matrix

[0228] Let the deterministic forecast mean field output by the decoder be μ1∈R. H×W×C Where H and W are the spatial grid sizes, and C is the number of predicted environmental factor categories. The covariance matrix Σ∈R of the covariance prediction head output is... (HW)×(HW) This is a high-dimensional matrix where the diagonal elements represent the variance (i.e., the magnitude of uncertainty) of the predicted values ​​of environmental factors at each spatial location, and the off-diagonal elements represent the correlation of prediction errors between different spatial locations or different factors. This matrix, together with the mean field, defines the multivariate Gaussian distribution of the atmospheric state at future times:

[0229] P(Y T+1 |X)=N(μ1,Σ).

[0230] (III) Efficient Modeling of Covariance Matrix

[0231] Directly outputting the complete (HW)×(HW) covariance matrix is ​​computationally and storage-infeasible. Therefore, this embodiment employs the same constraint modeling approach as Embodiment 5, combining block diagonalization with low-rank global correlation. For specific implementation details, please refer to the relevant description in Embodiment 5, which will not be repeated here.

[0232] (iv) Output and application of uncertain information

[0233] Through the above steps, this method ultimately outputs two types of results:

[0234] Deterministic forecast mean field μ1: the most likely concentration value of each factor at each location;

[0235] Covariance matrix Σ: describes the uncertainty of the forecast results and their spatial correlation.

[0236] Based on these two types of results, various derivative risk indicators can be further calculated:

[0237] Prediction confidence interval: For any spatial location and any environmental factor, calculate the prediction interval (e.g., 90% confidence interval) at a given confidence level.

[0238] Probability of exceeding the standard: For a given environmental quality standard threshold (e.g., daily average PM2.5 concentration of 75 μg / m³), 3 ), calculate the probability that the concentration at that location exceeds the threshold;

[0239] Spatial correlation analysis: Using the off-diagonal elements of the covariance matrix, the degree of correlation of forecast errors in different regions is analyzed to identify areas of systematic bias.

[0240] (v) Technical effects

[0241] The introduction of the uncertainty quantification step extends this method from deterministic forecasting to probabilistic forecasting, enabling the output to include not only "what is most likely to happen" but also "how probable it is." This extension provides crucial prior information for subsequent refined local site forecasting (see Example 6): the covariance matrix can serve as a priori input to the local model, allowing local probability estimation to perceive the confidence level of wide-area forecasts and achieving cross-scale transmission of uncertainty. Furthermore, in risk-sensitive scenarios such as environmental emergency decision-making and pollution early warning issuance, derived indicators such as the probability of exceeding standards have more direct guiding significance than a single concentration value.

[0242] It should be noted that this step is optional and can be enabled or disabled based on actual application needs. When uncertainty information is not required, only the deterministic mean field output by the decoder can be used; when risk quantification support is required, the covariance prediction head outputting the complete probability distribution should be enabled.

[0243] Furthermore, the covariance matrix is ​​modeled using a constraint method that combines block diagonalization and low-rank global correlation. The block diagonalization part is used to model the local correlation between different factors within each spatial location, while the low-rank global correlation part is used to model the long-range correlation across spatial locations.

[0244] This modeling approach significantly reduces computational complexity while ensuring the positive definiteness of the covariance matrix, making the estimation of high-dimensional covariance matrices feasible in engineering. Specifically, the block diagonal part is implemented by independently modeling a covariance matrix between factors for each spatial grid point; the local matrices of all grid points constitute a block diagonal matrix. The low-rank global correlation part is implemented by multiplying a learnable low-rank matrix by its transpose; its rank is much smaller than the total number of spatial locations, used to capture large-scale spatial correlations. The final covariance matrix is ​​expressed as the sum of the block diagonal matrix and the low-rank global correlation matrix, with a diagonal perturbation term added to ensure positive definiteness. Detailed implementation of this modeling approach can be found in the description of Example 5, and will not be repeated here.

[0245] Furthermore, the method also includes a vertical observation data fusion step: acquiring vertical observation data, and mapping the non-uniformly distributed vertical observation data onto a unified vertical grid through a learnable projection module, thereby realizing joint modeling of the ground observation field and the vertical observation field.

[0246] Vertical observation data (such as meteorological elements or atmospheric environmental factor concentrations at different pressure levels) are typically distributed across non-uniform altitude or pressure layers. Directly using them as input can lead to dimensional inconsistencies and misinterpretations of physical meaning. This step utilizes a learnable vertical field dynamic projection module to project the original non-uniform vertical layer features onto a pre-defined uniform vertical grid. The projection process can employ an attention-based mechanism or a learnable interpolation matrix. The learnable parameters are jointly optimized with the entire prediction model through end-to-end training, ensuring that the projected vertical features are most beneficial for the final atmospheric environmental factor forecast. The projected vertical features are then concatenated with the ground observation field parameters along the channel dimension to form a complete input feature containing vertical information, which is then used by the subsequent 3D spatiotemporal sensing encoder. Specific implementation details of this vertical observation data fusion step can be found in the description of the "Vertical Field Dynamic Projection Step" in Example 5, and are only briefly explained here.

[0247] After the above steps, this method ultimately outputs two types of results:

[0248] Deterministic forecast mean field: the most likely concentration values ​​of each factor at each location;

[0249] Covariance matrix: describes the uncertainty of forecast results and their spatial correlation.

[0250] Based on these two types of results, various derived risk indicators can be further calculated: prediction confidence interval, probability of exceeding the standard, etc. This uncertainty quantification result provides important prior information for subsequent local site-specific refined prediction (see Example 6)—the covariance matrix can be used as a priori input of uncertainty into the local model, enabling the local probability estimation to perceive the confidence level of the wide-area forecast and realize the cross-scale transmission of uncertainty.

[0251] The multimodal fusion method for predicting wide-area atmospheric environmental factors provided in this embodiment, through the collaborative design of a three-dimensional spatiotemporal sensing encoder, a multi-head joint decoder, and an optional uncertainty quantification module, has the following technical features:

[0252] (1) Computational efficiency of three-dimensional spatiotemporal modeling

[0253] Traditional numerical simulation methods (such as WRF-Chem and CMAQ) can generate continuous forecast fields, but they are computationally expensive, requiring several hours for a single high-resolution forecast, which limits their application in scenarios with high real-time requirements. This embodiment employs a 3D window attention mechanism based on a 3D visual Transformer, expanding the two-dimensional spatial window into a three-dimensional window that includes a time dimension. Self-attention computation is performed within the non-overlapping 3D window, simultaneously capturing the dynamic correlations of the atmosphere in horizontal space, vertical stratification, and time. Compared to models based on two-dimensional CNNs, this method can explicitly model physical correlations in the vertical direction (such as mass exchange between the boundary layer and the troposphere), enabling the capture of long-distance pollutant transport over thousands of kilometers. Furthermore, the windowed attention mechanism reduces computational complexity from the global quadratic level to the window level, making high-resolution wide-area forecasts feasible in engineering applications. The inference speed is also improved compared to traditional numerical models; for example, a 120-hour forecast can be completed in minutes.

[0254] (2) Compatibility and fusion of multi-source input data

[0255] This method can be used independently for forecasting using field spatial data such as meteorological reanalysis fields and atmospheric composition analysis fields, or it can input cross-modal spatiotemporal fusion features processed by Example 1. When using fusion features as input, sparse monitoring station observation data are injected into the field spatial data through end-to-end fusion, enabling wide-area forecast results to incorporate constraint information from local real-world observations and correct for systematic biases that may exist when relying solely on field data. This input design allows the method to adapt to application scenarios under different data conditions.

[0256] (3) Multi-factor collaborative prediction and parameter optimization

[0257] Through a multi-head joint output structure, the decoder shares the underlying feature layer and simultaneously learns the concentration forecast fields of various atmospheric environmental factors, utilizing the physical correlations between different factors, such as the correlation between PM2.5 and PM10, and the chemical coupling relationship between O3 and NO2. Compared with building a separate model for each factor, this method reduces the number of parameters by more than 50% and doubles the inference speed. Furthermore, due to the mutual complementarity of multi-task gradients, it exhibits certain generalization ability in sparse factors or regions. In addition, the efficient parameter fine-tuning strategy allows freezing the backbone network during deployment and fine-tuning only the output heads of specific factors, enabling the model to adapt to new regions or new pollution problems. The computational cost required for fine-tuning is reduced, for example, by more than 90%.

[0258] (4) Embedding of physical laws and consistency constraints

[0259] A physics-guided machine learning paradigm was introduced during training. Differentiable physical loss terms, such as mass conservation constraints and chemical kinetic constraints, enabled the model to learn to follow atmospheric physicochemical principles while fitting observational data. Experiments show that, with the introduction of physical constraints, the model's forecasts in unobserved areas exhibit more reasonable spatial distribution and temporal evolution patterns, suppressing physical inconsistencies that may occur in purely data-driven methods, such as anomalous pollutant accumulation in source-free areas and discrepancies between transport direction and wind field. Simultaneously, physical constraints act as a regularization mechanism, enhancing the model's generalization ability and robustness under extreme event scenarios.

[0260] (5) Uncertainty quantification capability

[0261] This method optionally integrates a covariance prediction head to output the covariance matrix of the forecast results, extending deterministic forecasts to probabilistic forecasts. By combining block diagonalization with low-rank global correlation in covariance modeling, it quantifies the prediction uncertainty of each factor at each location and the spatial correlation between different locations, while ensuring computational feasibility. This uncertainty information can be used to generate derived risk indicators such as exceedance probability and confidence intervals, supporting environmental emergency decision-making and pollution early warning issuance; it can also be used as prior information input into subsequent local refined prediction models (see Example 6) to achieve cross-scale transmission of uncertainty.

[0262] (6) Fusion of vertical observation data

[0263] To address the issue of uneven distribution of vertical observation data, this method designs a learnable vertical field dynamic projection module to map the non-uniform original vertical layer features onto a unified vertical grid. By optimizing the projection parameters through end-to-end training, the model can adaptively learn the physical relationships between different height layers, supplementing the 3D spatiotemporal modeling with vertical observation information. Experiments show that introducing vertical observation data fusion improves the model's simulation accuracy of the vertical distribution and transport processes of pollutants.

[0264] In summary, this embodiment provides background field support for wide-area atmospheric environmental factor forecasting for regional pollution prevention and control, emergency decision-making, and local refined forecasting through designs such as three-dimensional spatiotemporal sensing coding, multi-source input fusion, multi-factor collaborative prediction, physical constraint embedding, and uncertainty quantification.

[0265] Example 3

[0266] This embodiment relates to a dynamic spatial context encoding method, which belongs to the spatial feature extraction stage of the local model in the dual-expert fusion prediction architecture of the present invention. It is used to capture the spatial correlation between stations that changes in real time with wind direction and wind speed, and to provide high-quality spatial features for subsequent time series prediction (Embodiment 4).

[0267] In the field of spatiotemporal series prediction, the encoding of spatial context information is crucial to model performance. Currently, graph neural networks (GNNs) are widely used to model spatial relationships between sites. However, most GNN methods are based on static graph structures, meaning the edges and weights of the graph are predefined based on geographical distance or historical correlation. This static graph assumes that spatial relationships are fixed and cannot adapt to the objective laws governing the real-time changes in spatial influences in the atmospheric environment with conditions such as wind direction, leading to a decrease in prediction accuracy when dealing with sudden changes in meteorological conditions.

[0268] Based on this, this embodiment provides a dynamic spatial context coding method for air quality prediction, applicable to spatiotemporal prediction systems containing multiple monitoring stations. Specifically, it can be applied to scenarios requiring simultaneous consideration of spatial correlation and temporal dependence, such as air quality prediction, traffic flow prediction, and meteorological element forecasting. The following explanation uses air quality prediction as an example; however, those skilled in the art should understand that this method is equally applicable to other spatiotemporal prediction tasks.

[0269] like Figure 3 As shown, the method includes the following steps:

[0270] S300, Acquisition Steps: Acquire the multidimensional input features of all monitoring stations within the target area at the current discrete time step.

[0271] The monitoring stations are fixed or mobile observation points deployed within the target area to collect various types of environmental data, such as air quality monitoring stations, traffic flow detectors, and meteorological observation stations. In this embodiment, it is assumed that there are N monitoring stations in the target area, where N is an integer greater than or equal to 1.

[0272] The discrete time step refers to a time point obtained by sampling continuous time at fixed time intervals, such as every hour on the hour or every 15 minutes. In this embodiment, the time step is in hours, meaning each time step t corresponds to an hour on the hour. It should be understood that the time step size can be flexibly adjusted according to actual application requirements.

[0273] The multidimensional input features refer to the various observational data collected or acquired by each monitoring station at each time step t. For air quality prediction scenarios, multidimensional input features may include, but are not limited to:

[0274] Atmospheric environmental factors include: PM2.5 concentration, PM10 concentration, sulfur dioxide (SO2) concentration, nitrogen dioxide (NO2) concentration, ozone (O3) concentration, and carbon monoxide (CO) concentration.

[0275] Meteorological elements: such as temperature, humidity, air pressure, wind speed, wind direction, precipitation, etc.

[0276] Other supporting information includes timestamps (hour, day of the week, whether it is a holiday), emission source intensity in the area where the site is located (such as traffic flow, industrial emission index), etc.

[0277] The multidimensional input features also include the static attribute features of the monitoring station, such as the station's geographical location, surrounding environment type, and monitoring equipment attributes.

[0278] Specifically, the multidimensional input features also include location forecast information extracted from the wide-area background field. The wide-area background field is generated by the wide-area atmospheric environmental factor concentration forecasting method based on point-field data fusion as described in Embodiment 2 of this invention, and its output is a coarse-resolution atmospheric environmental factor concentration forecast field. For each monitoring station, the forecast values ​​of the grid where the station is located and its neighborhood are extracted from the coarse-resolution forecast field using bilinear interpolation or a local aggregation method based on an attention mechanism, serving as additional input features for the station at the current time step. This location forecast information characterizes the pollution transport background and meteorological forcing conditions in the station's region at a wide-area scale, providing macroscopic constraints for subsequent dynamic spatial context coding. The above data can be collected in real time by sensors at the monitoring station or obtained from external data sources such as meteorological departments and environmental monitoring centers. For missing data, conventional data preprocessing methods such as interpolation, padding, or using preceding time-series data can be employed.

[0279] At time t, the multidimensional input features of all N detection stations can be represented as a set X. t ={x 1t x 2t , ..., x jt , ..., x Nt}, where x jt Let x represent the multidimensional input feature vector of the j-th monitoring station at time step t, where j takes values ​​from 1 to N. Each x jt It is an F1-dimensional vector, where F1 is the dimension of the input features (i.e., the number of features), which is determined by the actual number of features selected.

[0280] This step provides the foundational data for subsequent dynamic spatial context coding, ensuring that observation information from all stations at each time step is fully input into the system.

[0281] It should be noted that although this embodiment uses an actual monitoring station as an example, the method is also applicable to spatial context encoding of any virtual location (i.e., sub-pixel) within the region. In this case, the monitoring station can be understood as any spatial location to be predicted, and its input features can be obtained through interpolation or extraction from a wide-area field.

[0282] S310, First transformation step: Perform a parameterized first transformation on the input features of each monitoring station to generate the source representation of the monitoring station at the current time step.

[0283] S320, Second transformation step: Perform a parameterized second transformation on the input features of each monitoring station to generate a content representation of the monitoring station at the current time step.

[0284] S330, Query generation step: Perform a parameterized third transformation on the input features of each monitoring station that serves as a target station to generate a query representation of the target station at the current time step.

[0285] In this invention, the first transformation step, the second transformation step, and the query generation step are feature transformation steps executed in parallel. The three share the same transformation rules and preprocessing requirements, and only differ in the transformation target, representation location, and execution object. The general transformation rules of the three are explained first, and then the differences in the implementation details of each step are described.

[0286] The general transformation rules for S310-S330 are as follows: The above three steps are feature transformation stages executed in parallel. They process the multi-dimensional input features of monitoring stations that have been collected and preprocessed (normalized, missing value imputation, etc.) in the acquisition step (S300). All three steps employ linear / nonlinear transformation methods learned autonomously during model training. The transformation parameters are model parameter matrices that can be updated through training iterations. There are no predefined fixed transformation rules, allowing for adaptive optimization based on the task objective and data features predicted from the spatiotemporal sequence. The transformation parameters are shared across all monitoring / target stations and all time steps, effectively controlling model complexity and improving generalization ability. Linear transformations are implemented using fully connected layers, achieved by adding bias terms through matrix multiplication of feature vectors and learned weight matrices. Nonlinear transformations are performed by applying activation functions such as ReLU, GELU, and Tanh after the linear transformation to achieve nonlinear mapping of features.

[0287] In this embodiment, all three preferably use the linear transformation ReLU nonlinear activation function to enhance the representation ability. When using a pure linear transformation, the activation function can be omitted. The specific choice is made by those skilled in the art based on the actual application requirements and model design.

[0288] If the multidimensional input features include static attribute features, all three are processed using the same fusion method: ① Joint encoding: After concatenating the static attribute features and dynamic observation features into a unified feature vector, a parameterized transformation is performed as a whole; ② Embedding fusion: First, the discrete static attribute features are mapped into a low-dimensional dense embedding vector through an independent embedding layer, and then the embedding vector is added or concatenated with the dynamic observation feature vector element by element before a parameterized transformation is performed.

[0289] When the multidimensional input features include location prediction information extracted from the wide-area background field, this information, together with dynamic observation features and static attribute features, constitutes a complete input feature vector. In the first transformation step, the second transformation step, and the query generation step, the location prediction information and the dynamic observation features are processed using the same parameterized transformation rules; that is, they are either directly concatenated and input into the transformation layer, or mapped through an independent feature encoding branch and then fused with other features. Specifically:

[0290] Direct concatenation and fusion: The location prediction information vector is concatenated with the dynamic observation feature vector and the static attribute feature vector to form the fused input feature, which is then input into the corresponding parameterized transformation layer to generate the representation;

[0291] Branch coding fusion: The dynamic observation features, static attribute features, and location forecast information are independently parameterized and transformed to generate corresponding intermediate representations. The three are then added or concatenated element by element and input into the subsequent transformation layer to achieve progressive fusion of multi-source information.

[0292] The parameters of the above fusion methods are all learned autonomously through model training, without the need for manual pre-setting of fusion weights. By introducing location prediction information of the wide-area background field, the source representation, content representation, and query representation not only include local observation information but also incorporate regional-scale pollution transmission background, enabling the subsequent dynamic spatial aggregation weight calculation to perceive wide-area environmental constraints.

[0293] Meanwhile, all three can be configured to execute multiple sets of independent parameterized transformations in parallel according to actual modeling needs. Multiple sets of corresponding transformations with different parameters and independent of each other are executed simultaneously on the same input feature. Each set of transformations generates a corresponding representation, which provides a basis for subsequent spatial interaction operations of multiple attention heads. The generated representation dimension is uniformly preset to the hyperparameter F2 (which can be adjusted according to the task complexity, such as 128 or 256), ensuring the dimensionality matching of subsequent inter-site interaction operations.

[0294] The first transformation step is implemented as follows:

[0295] Execution target: the multidimensional input features of each monitoring station within the target area at the current discrete time step. The feature transformation process of each monitoring station is independent of each other, ensuring the specificity of single-site feature processing and the parallelism of multi-site processing.

[0296] Characterization calculation: For monitoring station j at time step t, the input feature x jt Its source characterization k jt ∈R F2 Calculated using the following formula:

[0297] k j,t =Activation(W k ×xjt +b k );

[0298] Among them: W k ∈R F2 × F1 , is the learnable weight matrix; b k ∈R F2 , where is the learnable bias vector, and Activation(·) is the activation function.

[0299] Representation localization: Source representation k jt The aim is to capture the "queryable information" of monitoring station j at time step t, that is, the portion of the station's status information at the current moment that is available for other stations to query and focus on. In subsequent interactive computation steps, the source representations of all stations will jointly constitute a "source information database" for the target station's query representation to be retrieved and matched.

[0300] Output requirements: The source representation dimension generated by each monitoring station at the current time step should remain consistent to ensure the feasibility of subsequent inter-station interactive operations.

[0301] The second transformation step is implemented as follows:

[0302] Execution target: the multidimensional input features of each monitoring station within the target area at the current discrete time step, and the feature transformation processes of each monitoring station are independent of each other.

[0303] Parameter characteristics: The transformation parameter matrix in this step is completely independent of the first transformation step S110 and the query generation step S130, ensuring the uniqueness of the feature expression of the content representation.

[0304] Characterization calculation: For monitoring station j at time step t, the input feature x jt Its content represents v jt ∈R F2 Calculated using the following formula:

[0305] v j,t =Activation(W v ×x jt +b v );

[0306] Among them: W v ∈R F2 × F1 , is the learnable weight matrix; b v ∈R F2 , which is a learnable bias vector.

[0307] Representation localization: Content representation v jtThe aim is to capture the "information to be aggregated" of monitoring station j at time step t, that is, the core features of the current status information of the station that can be aggregated by the target station, which is the core feature carrier of subsequent spatial context encoding.

[0308] Output requirements: The content representation dimension generated by each monitoring station at the current time step should be consistent with the source representation dimension, and the content representation dimension of all stations should be unified to ensure the smooth execution of subsequent spatial interaction operations between stations.

[0309] The core difference between the first and second transformation steps lies in functional positioning and parameter independence:

[0310] The functions are different: the source representation generated in the first transformation step plays the role of "key," aiming to characterize the "information to be queried" at the current moment of the monitoring station, which is used to be retrieved and matched by the query representations of other stations, determining "where" the attention mechanism "focuses." The content representation generated in the second transformation step plays the role of "value," aiming to capture the "information to be aggregated" at the current moment of the monitoring station, serving as the core feature carrier for subsequent spatial context encoding, determining "what to take."

[0311] Parameter independence: The two steps use independent learnable weight matrices and bias vectors, without sharing parameters. This ensures that the model can learn different feature transformations for "retrievable representations" and "aggregable representations," thereby improving the richness and flexibility of the representation.

[0312] The subsequent usage differs: the source representations of all sites together constitute the "source information database," which is then used to calculate attention weights by performing a dot product with the query representation of the target site; while the content representations are weighted and summed according to the calculated weights to generate the final aggregated output.

[0313] In short, the two together constitute the "key-value pair" in the attention mechanism, and neither can be omitted, but their roles and parameters are completely independent.

[0314] The query generation process is implemented as follows:

[0315] Target site definition: refers to the monitoring site that performs spatial context encoding with respect to the current discrete time step; each monitoring site in the target area will be used as the target site at the corresponding time step to independently complete the generation of query representation, thereby realizing the modeling of spatial interaction relationships by enumerating all sites one by one.

[0316] Execution target: Multidimensional input features of all monitoring stations selected as target stations within the target area at the current discrete time step; the query generation process of each target station is independent of each other.

[0317] Parameter characteristics: The transformation parameter matrix in this step is completely independent of the parameter matrices in the first transformation step S110 and the second transformation step S120, ensuring that the feature representation of the query is unique and adapting to the dynamic calculation requirements of subsequent spatial aggregation weights.

[0318] Characterization calculation: For monitoring station j at time step t, the input feature x jt Its content represents q jt ∈R F2 Calculated using the following formula:

[0319] q j,t =Activation(W q ×x jt +b q );

[0320] Among them: W q ∈R F2 × F1 , is the learnable weight matrix; b q ∈R F2 , which is a learnable bias vector.

[0321] Representation localization: Query representation q jt The goal is to capture the query needs of monitoring site j at time step t, that is, the feature vectors used by the site at the current moment to retrieve source representations of other sites, which is the core basis for subsequent similarity calculation and spatial aggregation weight generation.

[0322] Output requirements: The query representation dimension generated for each target site at the current time step should be consistent with the source representation and content representation dimensions, and the query representation dimension of all target sites should be unified to ensure the smooth execution of similarity calculation and spatial aggregation weight generation in subsequent interactive operation steps.

[0323] S340, Interactive computation step: For each target site, interactive computation is performed using the query representation of the target site and the source representation of all monitoring sites to dynamically generate a set of spatial aggregation weights for the target site. Each of the spatial aggregation weights represents the aggregation weight of the target site for the information of a corresponding monitoring site at the current time step.

[0324] This step is the core computational step in the dynamic spatial context encoding method for mining time-varying spatial associations between stations. It follows the source representation output of the first transformation step (S310) and the query representation output of the query generation step (S330), providing accurate spatial attention weights for subsequent encoding generation steps. The execution targets are the query representation of each target station and the source representations of all monitored stations (including the target station itself) within the target area. The interactive computation processes of each target station are independent of each other, and the computation process does not require a predefined fixed spatial adjacency matrix. It can adaptively generate weights according to the actual correlation degree of station features, achieving accurate capture of dynamic spatial dependencies.

[0325] In this invention, the interactive operation is an attention-based operation. Through similarity calculation and normalization, the target site learns the spatial association strength with other monitoring sites autonomously. The weight value is positively correlated with the spatial influence of the corresponding monitoring site on the target site. The larger the value, the higher the contribution of the information of the monitoring site to the spatial context encoding of the target site.

[0326] The interactive operation step further includes two core sub-steps: similarity calculation and normalization processing. The specific execution process is as follows:

[0327] (1) Calculation of site association similarity: For any target site i, calculate its query representation q at the current time step t. it Source characterization k at the same time step t as all monitoring stations j within the target area jt Similarity is calculated enumerated one by one to obtain the set of similarity scores between the target site and all monitored sites. In this embodiment, the scaled dot product method is preferably used to calculate the similarity, which can effectively alleviate the problem of excessively large inner product values ​​under high-dimensional representation. The specific calculation formula is as follows:

[0328] sim(q it k jt ) = (q jt ·k jt ) / F 1 / 2 Where, q jt ·k jt F is the vector inner product of the query representation and the source representation, where F is the representation dimension (consistent with the preset dimensions of the source and query representations). Those skilled in the art can replace it with other conventional similarity calculation methods such as cosine similarity or dot product similarity according to actual modeling needs.

[0329] Normalization: All similarity scores obtained above are normalized, mapping them to the 0-1 interval, with the sum of all scores being 1. This ensures the results satisfy probability distribution characteristics, facilitating subsequent weighted aggregation operations. In this embodiment, the Softmax function is preferably used for normalization. The specific calculation formula is as follows:

[0330] .

[0331] Where, α ijt The normalized attention weights, i.e., the spatial aggregation weights, are represented by r, which is the enumeration number of the monitoring station. After this step, each target station i will generate a set of N-dimensional spatial aggregation weights α. ijt =[α i1t , ..., α ijt , ..., α iNt ].

[0332] To address the modeling requirements of multi-head attention, the interactive computation steps are further configured as follows:

[0333] Based on the multiple sets of source representations generated in the first transformation step and the multiple sets of query representations generated in the query generation step, each set is matched one by one. For each set of matched query representations, an independent attention-weighted aggregation process (i.e., the aforementioned similarity calculation and normalization) is performed in parallel with the source representations. Each set of operations independently generates a set of spatial aggregation weights, ultimately generating multiple different sets of spatial aggregation weights for each target site. These multiple sets of spatial aggregation weights capture the spatial interaction relationships between the target site and other monitoring sites from different dimensions, effectively mining complex, multi-dimensional spatial correlation features between sites. This provides dedicated weight support for the subsequent generation of multiple sets of intermediate context codes, improving the model's ability to express complex spatiotemporal data.

[0334] After the above interactive processing, each target site generates one or more sets of spatial aggregation weights. The dimension of each set of spatial aggregation weights is consistent with the total number of monitoring sites in the target area, and the dimension of the spatial aggregation weights generated by all target sites at the same time step is consistent. In the multi-set mode, the dimension of each set of weights is consistent with the single-set mode, and each set of weights is independent of each other, ensuring the smooth execution of weighted aggregation operations in the subsequent encoding generation steps.

[0335] S350, Encoding generation step: For each target site, according to the spatial aggregation weight corresponding to the target site, the content representation of all monitored sites is weighted and aggregated to generate the spatial context encoding of the target site at the current time step.

[0336] This step is the core link in the dynamic spatial context coding method to realize spatial information fusion. It follows the spatial aggregation weight output of the interactive operation step (S340) and the content representation output of the second transformation step (S320). Through weighted aggregation, it fuses the spatial association information of the target station with all other monitoring stations and the core features of each station, and finally generates a context code that can represent the spatial interaction relationship of the target station at the current time step. The execution objects are the spatial aggregation weight corresponding to each target station and the content representation of all monitoring stations (including the target station itself) in the target area. The coding generation process for each target station is independent of each other, ensuring the uniqueness and specificity of the spatial context code of a single station.

[0337] The core execution logic of this step is weighted summation aggregation. That is, using the spatial aggregation weights generated in the interactive computation steps, the content representations of all monitoring sites are linearly weighted, ensuring that the contribution of the content representations is positively correlated with the corresponding spatial aggregation weights. The larger the spatial aggregation weight, the higher the proportion of the content representation of the corresponding monitoring site in the spatial context encoding of the target site. This accurately characterizes the spatial influence of different monitoring sites on the target site. For any target site i, the spatial context encoding c at the current time step t... it Calculated using the following formula:

[0338] .

[0339] c it The dimensions of the feature should be consistent with the dimensions of the content representation to ensure the uniformity of the feature dimensions.

[0340] To address the modeling requirements of multi-head attention, this step is further configured as an execution mode of multiple sets of parallel weighted aggregation. Specifically, based on the multiple sets of spatial aggregation weights generated for each target site in the interactive operation step, and the multiple sets of content representations generated for all monitoring sites in the second transformation step, these sets are matched one by one, with each set of matching being a set of spatial aggregation weights corresponding to a set of content representations. For each set of matched weights and representations, the above weighted sum aggregation operation is performed independently, and each set of operations generates a corresponding intermediate context code, ultimately generating multiple different intermediate context codes for each target site.

[0341] In the process of multi-group parallel weighted aggregation, the operations of each group are independent of each other. The weighted summation of each group is completed based on the exclusive spatial aggregation weight and content representation. Different intermediate context codes can capture the spatial fusion features of the target site and other monitoring sites from different dimensions, and fully explore the complex and multi-dimensional spatial interaction relationships between sites. Moreover, the dimensions of the intermediate context codes generated by each group are consistent with the dimensions of the spatial context codes in the single-group mode. All intermediate context codes are unified, providing a standardized and adapted feature foundation for feature integration in the subsequent multi-head fusion steps.

[0342] After the above weighted aggregation process, each target site will generate a spatial context code (single-head attention mode) or multiple different intermediate context codes (multi-head attention mode). The coding dimensions generated by all target sites at the same time step are consistent, and the coding features fully integrate the spatial association information between the target site and all monitoring sites in the region. This can accurately characterize the spatial environment and interaction relationship of the target site at the current time step, providing high-quality spatial coding features for the subsequent time series prediction module.

[0343] Furthermore, the method also includes a post-processing step:

[0344] S360, the spatial context code is input into one or more nonlinear transform layers for processing to generate the final spatial coded features.

[0345] This step is the feature enhancement and dimension adaptation stage of spatial context coding. It follows the output of the coding generation step and improves the expressive power of the features through nonlinear transformation, so that the spatial coding features are more suitable for the input requirements of the subsequent spatiotemporal sequence prediction model. The execution targets are the spatial context codes generated by the coding generation step for each target site (single-head mode) or the final spatial context codes obtained by the multi-head fusion step (multi-head mode), and the post-processing processes for each target site are independent of each other.

[0346] The nonlinear transformation layer can be a combination of fully connected layers and nonlinear activation functions, or it can be a multilayer perceptron (MLP) structure. In this embodiment, a combination of "fully connected layer + GELU activation function + Dropout layer" is preferred, which ensures the nonlinear mapping capability of features and mitigates the risk of model overfitting. Specifically, for the spatial context encoding c of target station i at the current time step t... it Its final spatial coding feature E it Calculated using the following formula:

[0347] E it =Dropout(GELU(W e ×C it +b e )).

[0348] Among them, W e ∈R FE×F2 , is a learnable transformation weight matrix, and FE is the dimension of the final spatial encoded features (which can be preset according to the needs of subsequent prediction tasks, such as matching the dimension of the prediction target).

[0349] b e ∈R FE , is a learnable bias vector;

[0350] GELU(·) is a non-linear activation function that is better suited to the smooth mapping of high-dimensional features than ReLU. It can also be replaced by ReLU, Swish, etc., as needed.

[0351] Dropout(·) is a random deactivation layer with a deactivation probability that can be set to 0.1~0.3 (0.2 is preferred in this embodiment) to improve the generalization ability of the model.

[0352] If multiple nonlinear transformation layers are configured (such as two fully connected layers), a layer-by-layer progressive transformation method is adopted: c1 it =GELU(W e1 ×c it +b e1 ); E it =Dropout(GELU(W e2 ×c1 it +b e2 )).

[0353] Among them, W e1 ∈R FM×F2 W e2 ∈R F3×FM FM represents the weight matrix for different layers, and FM is the intermediate transformation dimension (which can be set to the middle value between F2 and F3, such as 256) to ensure the smoothness of feature transformation.

[0354] After the post-processing step, the final spatial coding feature dimension of each target site is unified, and the feature incorporates high-order spatial interaction information after nonlinear transformation, which can effectively improve the ability of subsequent spatiotemporal sequence prediction models to capture complex spatial patterns. At the same time, the dimension of the final spatial coding feature can be flexibly adapted to the input requirements of different prediction tasks (such as single-point concentration prediction and regional trend prediction), enhancing the versatility of the method.

[0355] Furthermore, the post-processing step also includes a residual connection: adding the spatial context code or the derived features of the spatial context code to the output of the nonlinear transform layer element by element to generate the final spatial code features.

[0356] To alleviate the gradient vanishing and network degradation problems that may occur during deep network training and to improve the stability and representation ability of feature transformation, this embodiment provides a post-processing implementation method that includes residual connections; at the same time, layer normalization can be introduced to further optimize the training process and ensure the model convergence efficiency and stability.

[0357] The core of residual connections is to construct a parallel structure of identity mapping branches and nonlinear transformation branches, allowing gradients to propagate directly back through the identity mapping branches, thus overcoming the training bottleneck of deep networks. In practice, it is necessary to ensure that the input and output dimensions of the residual branches are strictly consistent to avoid computational errors caused by dimension mismatch.

[0358] Spatial context encoding c for target site i at the current time step t it (Dimension F2), the basic implementation of residual connections is: s it =c it +f(c it ).in:

[0359] f(·) is the transformation function of the aforementioned nonlinear transformation layer (such as a multilayer perceptron MLP), and its output dimension needs to be the same as c. it Maintain consistency to ensure that element-wise addition operations are feasible;

[0360] s it These are the initial spatial coding features after residual connection.

[0361] To further enhance feature representation, a lightweight nonlinear activation function (such as GELU or ReLU) can be added after the residual summation, in the form: s it =Activation(c it +f(c it Among them, the GELU function is preferred for Activation(·), which balances smoothness and non-linearity. Those skilled in the art can choose whether to add this activation function according to the model depth (shallow networks can omit it, while deep networks are advised to retain it).

[0362] The MLP with nonlinear transform branches needs to ensure that the input and output dimensions are consistent. In this embodiment, a structure of two fully connected layers + activation function + Dropout is preferred, specifically: f(c it =Dropout(GELU(W2×GELU(W1×c)) it +b1)+b2))where:

[0363] W1∈R Fh×F2 W2∈R F2×Fh Fh is the learnable weight matrix, and Fh is the dimension of the intermediate hidden layer (which can be set to 2×F2, such as Fh=256 when F2=128).

[0364] b1∈R Fh b2∈R F2 , is a learnable bias vector;

[0365] The deactivation probability of Dropout(·) is set to 0.2 to alleviate overfitting; at this time, the complete formula for residual connection can be simplified to: s it =c it +MLP(c it ), where MLP(·) is the transformation of the two-layer perceptron mentioned above, ensuring that the output dimension is the same as c. it Exact match.

[0366] To further stabilize the training process and accelerate model convergence, this embodiment introduces layer normalization at the input of the nonlinear transformation layer to standardize the feature distribution of spatial context encoding. Layer normalization independently calculates the mean and variance for the feature dimension of each sample, adapting to the sample-by-sample distribution characteristics of site-level features, and outperforming batch normalization in the case of small batches of data.

[0367] The residual connection implementation after introducing layer normalization is as follows: s it =c it +MLP(LayerNorm(c it The specific calculation method for LayerNorm(·) is as follows: LayerNorm(c it ) = (c it -μ it ) / [(σ it ) 2 +ε] 1 / 2 ×γ+β.

[0368] μ it For c it The mean of the feature dimensions, (σ it ) 2 Variance of feature dimension; ε=10 -5 To prevent the denominator from being a minimum value of 0; γ and β are learnable scaling and offset coefficients to ensure the flexibility of the normalized features.

[0369] After the above residual connection (optional layer normalization) processing, s can be transformed without additional transformation. it As the final spatial encoding feature; if it is necessary to adapt to the dimensionality requirements of the subsequent temporal prediction module, a dimensionality-adaptive linear layer (without an activation function) can be added after the residual connection:

[0370] E it =W out ×s it +b out .

[0371] Among them W out ∈R Fout×F2 bout ∈R Fout Fout is the input dimension required by subsequent prediction modules (e.g., 64, 128).

[0372] The final spatial coding features of all monitoring stations at all historical time steps can be integrated into a tensor with shape (N, Thistory, Fs):

[0373] Thistory is the length of a historical time step (e.g., the past 24 hours).

[0374] Fs is the dimension of the spatial encoding feature (i.e., F2 or Fout); this tensor fully preserves the dynamic spatial context information of each station at the historical time step, and serves as the core input of the subsequent time series prediction module to achieve spatiotemporal sequence prediction for future time steps (such as the next 12 / 24 hours).

[0375] Furthermore, the method also includes:

[0376] S370, Multi-head fusion step: The multiple intermediate context codes are concatenated or linearly transformed to obtain the final spatial context code.

[0377] This step is the feature integration stage of spatial context encoding under the multi-head attention mechanism. It takes multiple intermediate context codes output from the encoding generation step and integrates multi-dimensional spatial interaction features into a single final spatial context code through a unified fusion method. This not only preserves the multi-dimensional spatial correlation information mined by multi-head attention but also ensures the uniformity of feature dimensions, adapting to the input requirements of subsequent post-processing steps. The execution target is all intermediate context codes generated by the encoding generation step for each target site (the number is consistent with the number of attention heads / parallel transformation groups, denoted as m), and the multi-head fusion process for each target site is independent of each other.

[0378] I. Splicing and fusion method (dimensional expansion type)

[0379] Concatenation fusion preserves all feature information of the head by directly concatenating the feature dimensions of multiple sets of intermediate context codes, making it suitable for scenarios requiring maximum feature richness. Specifically, for m intermediate context codes (each with an F2 dimension) generated by target site i at the current time step t, the concatenation fusion calculation method is as follows: c final-it =Concat(c 1 it , ..., c 1 im )in:

[0380] Concat(·) is a dimension concatenation operation that is performed along the feature dimension (the last dimension).

[0381] After splicing cfinal-it The dimension is m×F2 (e.g., when m=8 and F2=128, the dimension after concatenation is 1024); to adapt to the dimension requirements of subsequent steps, a linear transformation layer needs to be added after concatenation to map the dimension back to F2, in the form: c final-it =W concat ×Concat(c 1 it , ..., c 1 im )+b concat Among them, W concat ∈R F2×(m×F2) b concat ∈R F2 These are learnable linear transformation parameters that ensure the fused dimension is consistent with the spatial context encoding in the single-head mode.

[0382] II. Linear Transformation Fusion Method (Dimension Preservation)

[0383] Linear transformation fusion integrates multiple sets of intermediate context codes through weighted summation, directly maintaining the feature dimension as F2 without requiring additional dimension adaptation. This is suitable for scenarios requiring control over feature dimensions and reduced computational complexity. The core logic assigns learnable weights to each attention head, achieving adaptive weighted fusion of multiple feature sets. The specific calculation method is: c final-it =(W h ×c 1 it )+……+(W h ×c h it )+……+(W h ×c m it ))+b liner .in:

[0384] h is the enumeration number of the attention head (1≤h≤m);

[0385] W h ∈R F2×F2 Let h be the learnable weight matrix for the h-th head. The weights for each head are learned independently to match the importance of the spatial correlation features captured by that head.

[0386] b linear ∈R F2 This is the global bias vector;

[0387] The summation result is directly used as the final spatial context encoding, with the dimension maintained at F2, which is fully compatible with the single-head mode.

[0388] A simplified implementation of linear transformation fusion is as follows: all intermediate context codes are stacked into a tensor c. stack-it ∈Rm ×F2 Blending is achieved through a shared linear projection layer: c final-it =W proj ×Flatten(c stack-it )+b proj Where Flatten(·) is the flattening operation (flattens an m×F2 tensor into an m×F2 vector), W proj ∈R F2×(m×F2) b proj ∈R F2 As learnable parameters, this form is easier to implement in engineering.

[0389] The advantage of stitching and fusion is that it fully preserves all spatial feature information discovered by multiple sources, with high feature richness, making it suitable for scenarios with a large number of monitoring stations and complex spatial relationships (such as city-level air quality prediction). The disadvantage is that the dimensions are temporarily expanded after stitching, requiring additional linear transformation and slightly higher computational cost.

[0390] The advantages of linear transformation fusion are that the dimensionality remains consistent, the computational efficiency is high, and the number of model parameters is small, making it suitable for lightweight deployments or scenarios with a small number of monitoring stations (such as county-level meteorological forecasting). The disadvantage is that some fine-grained features may be lost through weighted summation. Those skilled in the art can choose the corresponding fusion method according to the complexity of the actual task, computing resources, and prediction accuracy requirements. Both methods can effectively integrate multiple sets of intermediate context codes, and the fusion parameters are all learned autonomously through model training without the need for manual preset.

[0391] After the multi-head fusion process, the final spatial context code generated for each target site must meet the following requirements:

[0392] Dimensional uniformity: The final spatial context encoding dimension of all target sites at the same time step is F2, which is consistent with the spatial context encoding dimension in single-head mode;

[0393] Feature integrity: The fused encoding must fully preserve the spatial interaction relationships of different dimensions captured by multiple sets of intermediate context encodings, so as to provide a high-quality spatial feature foundation for subsequent post-processing steps;

[0394] Computational independence: The multi-head fusion process of each target site is independent of each other, avoiding cross-site feature interference and ensuring the specificity of single-site spatial context encoding.

[0395] Furthermore, the method also includes:

[0396] S380, Future Information Fusion Step: Obtain known or predicted external driving features of the target region in future time steps; when performing spatial context encoding on the current time step, use the external driving features of the future time steps as auxiliary information, and use them together with the input features of the current time step to generate the query representation, source representation, or content representation.

[0397] This step is an extended optimization of spatiotemporal feature fusion. Its core purpose is to break the problem of spatiotemporal information fragmentation caused by "relying solely on historical / current information for spatial encoding". By introducing external driving features from future time steps, spatial context encoding can capture the potential impact of future environmental changes on the spatial association of the current site in advance, thereby improving the adaptability of encoded features to future time series prediction. The execution timing is parallel to the acquisition step (S300), and feature fusion is completed before the generation of source representation, content representation, and query representation. The fusion process is uniformly executed for all monitoring sites / target sites.

[0398] The external driving features of the future time step are environmental features that are strongly correlated with the spatiotemporal sequence prediction task and can be obtained or predicted in advance, specifically including but not limited to:

[0399] Known deterministic features: artificial planning features such as statutory holiday markers, traffic control plans, and industrial emission scheduling in future time steps;

[0400] Predictive features include meteorological forecast data (wind speed, wind direction, temperature, precipitation) for future time steps, large-scale circulation characteristics of air quality forecasts, and traffic flow predictions. These features have prediction biases and require error calibration before use. The scope of the future time step should match the time span of subsequent prediction tasks (e.g., if predicting the next 24 hours, obtain the external driving features for the next 24 hours), and the time granularity should remain consistent with the current discrete time step (e.g., hourly).

[0401] For known deterministic features, they are directly collected from business systems / public data sources, and after completing discrete feature encoding (such as one-hot encoding and embedded encoding) and continuous feature normalization, a feature vector with unified dimensions is formed.

[0402] Furthermore, the external driving features of the future time step may also include covariance matrix information generated by the wide-area probabilistic forecasting method described in Embodiment 5 of the present invention. This covariance matrix information characterizes the prediction uncertainty of the wide-area model at different spatial locations. After compression encoding, it is fused with the query representation or source representation in this step, enabling the dynamic spatial context encoding to perceive the confidence level of the wide-area forecast and providing prior information for subsequent uncertainty estimation.

[0403] Future externally driven features need to be fused with the input features at the current time step before being input into the first / second transformation step and query generation step to generate the corresponding representation. The fusion needs to ensure dimensionality matching, temporal alignment, and feature complementarity. This embodiment provides two core fusion methods, which can be selected according to the feature type:

[0404] 1. Temporal aggregation and fusion (adapting to features of multiple future time steps)

[0405] First, perform temporal aggregation on the external driving features of future time steps to extract core trend information, and then fuse it with the input features of the current time step:

[0406] D agg =MeanPooling(D future ) or Dagg=MaxPooling(D future ).

[0407] Among them, D agg For the aggregated future-driving features, MeanPooling / MaxPooling uses mean / maximum pooling to aggregate the future features of TF×Fd into an aggregated feature vector of 1×Fd; then the aggregated features are combined with the input features x at the current time step. it splicing: x m it =Concat(x it D agg Where, x m it The input feature vector after fusion has a dimension of Fc+Fd (Fc is the dimension of the current input feature).

[0408] 2. Time-step alignment and fusion (adapting to single future time-step features)

[0409] If only the external driving features of a single future time step aligned with the current time step are obtained (e.g., if the current time step is t, obtain the external driving features at time t+1), the future feature vector d can be directly used. future,t Add or concatenate element-wise with the current input features:

[0410] splicing and blending: x it =Concat(x) it ,d future,t )

[0411] Element-wise addition: First, map the future feature dimensions to Fc using a linear transformation, then add them together: x m it =x it +W d ×d future,t +b d Among them, W d ∈RFc×Fd b d ∈R Fc These are learnable linear transformation parameters.

[0412] The fused input features x m it The original input features are replaced by the corresponding transformation steps to generate representations. The core rule is:

[0413] Unified transformation rules: The first transformation step, the second transformation step, and the query generation step all adopt the same parameterized transformation logic as the original steps, only replacing the input features with the fused x. m it The dimension of the transformation parameter matrix is ​​adapted to F2×(Fc+Fd) (splicing and blending) or kept at F2×Fc (element-by-element addition and blending).

[0414] Optional partial fusion: Based on the needs of the prediction task, only future external driving features can be incorporated into a certain type of representation (e.g., only query representations can be incorporated to capture the target site's query needs for the future environment). Specifically:

[0415] Source representation / content representation: still based on the original input features x it generate;

[0416] Query representation: based on the fused input features x m it This method allows for precise control over the scope of future information, avoiding overfitting caused by excessive introduction of future information.

[0417] The technical rationale and effectiveness of this step are as follows:

[0418] Temporal complementarity: Future external driving characteristics (such as wind direction forecast) directly affect the spatial correlation between monitoring stations (such as the direction of pollutant transport). Incorporating spatial coding in advance allows the model to capture the potential impact of future environmental changes on the current spatial correlation, solving the problem that spatial coding in existing technologies only reflects the current state and is disconnected from future prediction targets.

[0419] Feature adaptability: The fusion method takes into account both multi-time step aggregation and single-time step alignment, adapting to different types of future externally driven features, and the fusion parameters are learned autonomously through training, without the need for manual preset of fusion weights;

[0420] Risk control: Supports local fusion mode, which can avoid the decrease in spatial coding accuracy due to the error propagation of future predicted features, thus ensuring the robustness of the model.

[0421] This invention constructs a site-level dynamic spatial context coding system based on a self-attention mechanism, systematically designed from multiple levels including feature processing, spatial association modeling, multi-head parallel coding, residual stable training, and future information fusion. Compared with traditional fixed-topology graph neural networks, simple spatial weighting methods, and general attention coding methods, this method has the following technical features:

[0422] 1. Dynamic adaptive spatial dependency modeling, without the need for predefined adjacency matrices.

[0423] By performing attention-based interaction operations between query representations and source representations, the spatial correlation strength between sites can be learned in real time. Without relying on adjacency matrices constructed artificially based on geographical distance, road network structure, etc., it can capture time-varying spatial dependencies in spatiotemporal data that change with wind field, emissions, and meteorological conditions.

[0424] 2. Three-branch independent parameterized transformation to achieve decoupling of query-source-content features.

[0425] This method employs three independent learnable branches: a first transformation, a second transformation, and query generation. These branches generate source representations, content representations, and query representations, respectively, achieving functional separation of "what to focus on, what to aggregate, and what to retrieve." Compared to single transformations or simple concatenation methods, this approach differs in feature representation capability, spatial matching accuracy, and model generalization.

[0426] 3. Supports multi-head attention parallel encoding, mining spatial structure from multiple dimensions.

[0427] It supports multiple independent transformations, multi-head attention parallel computation and multi-head fusion, and can capture the spatial relationship between sites from different dimensions and scales at the same time. Compared with single attention heads, it can provide richer spatial context features and higher robustness.

[0428] 4. Introduce residual connections and layer normalization to improve the training stability of deep models.

[0429] In post-processing, the use of residual structure combined with layer normalization can alleviate the problems of gradient vanishing, network degradation and training instability in deep networks to a certain extent, enabling the model to maintain convergence speed and training stability in deep structures, while enhancing the nonlinear expressive power of high-order features.

[0430] 5. Supports joint encoding of static attributes and dynamic observation features to enhance site heterogeneity modeling.

[0431] By jointly encoding or embedding static attributes such as the geographical location, equipment type, and surrounding environment of monitoring stations with dynamic observation characteristics such as real-time weather and concentration, the spatial coding can simultaneously include spatial structure information, inherent station attributes, and real-time status information, thus adapting to heterogeneous monitoring networks.

[0432] 6. Introduce future external driving information fusion to enhance predictive guidance.

[0433] By incorporating known or predictable external driving features such as future weather, future emissions, and future planning into the spatial coding stage, the spatial context coding at the current moment can reflect future environmental change trends, which is different from methods that only use historical information for coding.

[0434] 7. The entire process is modularized and step-by-step, making it highly versatile and easy to deploy in engineering projects.

[0435] The spatial coding process is broken down into standardized steps such as acquisition, transformation, interaction, encoding, post-processing, and fusion. The structure is clear, the reproducibility is high, and it can be embedded as a general spatial encoder into various spatiotemporal sequence prediction models (such as air quality prediction, traffic flow prediction, and weather prediction), with a wide range of applications.

[0436] In this embodiment, the spatial context encoding is used for probabilistic downscaling and uncertainty estimation of atmospheric environmental factor concentrations at the station level. It is input into the parameterized network as an enhanced feature of historical observation data from surrounding monitoring stations, and the parameterized network outputs the probability distribution parameters for any predicted location.

[0437] Furthermore, the method provided in this embodiment can also be applied to probabilistic downscaling scenarios:

[0438] The spatial context codes generated by each monitoring station at multiple historical time steps are combined in chronological order to form a spatial context code feature sequence.

[0439] The spatial context encoding feature sequence is used as an enhanced feature of historical observation data from surrounding monitoring stations, and is input into the parameterized network along with wide-area background field data and static geographic attribute data.

[0440] Given the geographic coordinates of any predicted location within the target area, the parameterized network outputs a set of parameters representing the probability distribution of atmospheric environmental factor concentrations at the predicted location in future time steps.

[0441] Based on the probability distribution parameter set, a complete probability density function for the concentration of the atmospheric environmental factors is constructed.

[0442] For a detailed description of the above-mentioned probabilistic downscaling, please refer to Embodiment Six of the present invention.

[0443] The dynamic spatial context encoding generated in this embodiment can not only be directly input into the time series prediction module for site prediction, but also serve as an enhanced feature for historical site observation data, applied to the site-level probabilistic downscaling method described in Embodiment Six. By inputting the spatial context encoding feature sequence containing dynamic spatiotemporal interaction information between sites, along with wide-area background field data and static geographic attribute data, into a parameterized network, the probabilistic prediction at any location can fully benefit from a deep understanding of the spatiotemporal relationships of surrounding sites, thereby significantly improving prediction accuracy and the rationality of uncertainty estimation.

[0444] Based on the same inventive concept, embodiments of the present invention provide a dynamic spatial context coding apparatus, comprising:

[0445] The input interface is used to obtain the multidimensional input features of all monitoring stations within the target area at the current discrete time step;

[0446] The processing unit is configured to execute the aforementioned method to generate a corresponding spatial context code for each monitoring station that serves as a target station;

[0447] The output interface is used to output the spatial context encoding sequence of all monitoring stations at the current time step.

[0448] The input interface serves as the interaction point between the device and external data sources. Its core function corresponds to the feature acquisition step (S300) and the future information fusion step (S380) of the aforementioned method. On the one hand, the input interface collects multi-dimensional input features (including dynamic observation features and static attribute features) of all monitoring stations within the target area at the current discrete time step, and completes preprocessing such as feature normalization, missing value filling, and discrete feature embedding, outputting a feature vector with unified dimensions to provide compliant input for subsequent transformation steps. On the other hand, the input interface can be extended to access known / predicted external driving features (such as weather forecasts and traffic scheduling data) of the target area at future time steps, and completes feature calibration and format conversion according to the preprocessing rules defined in the method, providing data support for future information fusion.

[0449] The processing unit is the core computing module of the device. It integrates multiple sets of parameterized transformation submodules, interactive operation submodules, encoding generation submodules, and post-processing submodules. It is configured to completely execute the first transformation step (S310), the second transformation step (S320), the query generation step (S330), the interactive operation step (S340), the encoding generation step (S350), the post-processing step (including residual connection / layer normalization), and the multi-head fusion step (S370) of the aforementioned method.

[0450] Transformation submodule: Performs parameterized first / second / third transformations on the feature vectors output from the input interface to generate source representation, content representation, and query representation respectively. It supports multiple sets of parallel transformation modes to meet the modeling requirements of multi-head attention.

[0451] Interactive computation submodule: Based on the attention mechanism, it calculates the similarity score between the target site query representation and the source representation of all monitored sites, and generates dynamic spatial aggregation weights after normalization to realize adaptive mining of time-varying spatial associations between sites;

[0452] Encoding generation submodule: Utilizes spatial aggregation weights to perform weighted aggregation on content representations, generating intermediate context encoding (multi-head mode) or spatial context encoding (single-head mode) for the target site.

[0453] Post-processing submodule: Performs nonlinear transformations, residual connections, layer normalization, and other operations on the encoding results to alleviate the degradation problem of deep networks and output high-quality final spatial coding features;

[0454] Multi-head fusion submodule (optional): Performs concatenation / linear transformation fusion on multiple sets of intermediate context codes, outputting a final spatial context code with unified dimensions to meet the input requirements of subsequent time series prediction modules. All computational logic of the processing unit strictly follows the parameter definitions of the aforementioned methods (such as representation dimension F2, weight matrix dimension, activation function selection), and the computational processes of each submodule are independent of each other, ensuring the specificity of single-site encoding and the parallelism of multi-site processing.

[0455] The output interface, as the result output terminal of the device, corresponds to the result output stage of the encoding generation step and the post-processing step in the aforementioned method: the output interface integrates the spatial context codes (or the final spatial coding features after post-processing) of all monitoring stations generated by the processing unit at the current time step into an ordered tensor sequence, with a tensor shape of (N, Fs) (N is the total number of monitoring stations, and Fs is the dimension of the coding features); if it covers encoding of multiple historical time steps, the output interface can be expanded to have an output shape of (N, T) history The spatiotemporal coding tensor of Fs is directly used as the core input of the subsequent spatiotemporal sequence prediction module, achieving seamless integration with the prediction task.

[0456] This device, through hardware-level parallel computing optimization (such as GPU acceleration) and software-level parameter sharing mechanism design, not only ensures the complete implementation of the core innovations of the aforementioned methods (dynamic spatial aggregation weights, residual connections, future information fusion, etc.), but also takes into account coding efficiency and engineering practicality. It can be widely adapted to the site-level spatial context coding needs of spatiotemporal sequence prediction scenarios such as air quality, meteorology, and traffic flow.

[0457] Example 4

[0458] This embodiment relates to a local site hybrid expert temporal prediction method, which belongs to the temporal prediction link of the local model in the dual expert fusion prediction architecture of the present invention. It is used to perform temporal modeling on spatially encoded features, output site-level refined prediction results, and can be further used for probabilistic downscaling (Embodiment Six).

[0459] Temporal forecasting is a core component of air quality forecasting. Existing deep learning prediction models typically tightly couple spatial encoding and temporal forecasting modules into a single, integrated model. This tightly coupled design has several drawbacks: First, when trying new temporal forecasting algorithms (such as replacing LSTM with Transformer), the entire model often needs to be reconstructed, leading to high iteration costs and poor flexibility. Second, pollution temporal patterns at different monitoring stations may vary significantly (e.g., traffic stations are significantly affected by morning and evening rush hours, while background stations exhibit stable diurnal variations). A single model struggles to adapt to the local characteristics of all stations simultaneously, limiting its generalization ability and resulting in unstable prediction accuracy in monitoring networks with diverse station types.

[0460] To address the aforementioned technical issues, this embodiment provides a hybrid expert time series prediction method applied to local site prediction. The core concept of this method lies in: achieving the "pluggable" characteristic of the time series prediction module through standardized interface design, enabling the time series prediction algorithm to be flexibly replaced independently of the spatial coding module; simultaneously introducing a Mixture of Experts (MoE) architecture, utilizing a gating network to dynamically allocate expert weights based on site characteristics, allowing different site types or different pollution processes to be modeled by specialized expert sub-units, thereby achieving adaptive prediction of diverse time series patterns.

[0461] like Figure 4 As shown, the method includes the following steps:

[0462] S400, obtain the spatial context coding feature sequence that characterizes the spatial interaction of each monitoring station.

[0463] The spatial context-coded feature sequence is a temporal feature representation that integrates dynamic spatial correlation information between sites. Its core function is to incorporate wide-area transmission background and local spatial interaction information into subsequent temporal prediction stages. This feature sequence is generated by a dynamic spatial context coding method, specifically:

[0464] The dynamic spatial context coding method independently processes the multidimensional input features of all monitoring stations at each discrete time step. Through parameterized transformation and attention mechanisms, it generates a context code for each monitoring station that can characterize its current spatial environment. The above process is repeated for multiple historical time steps, and the context codes generated at each time step are combined in chronological order to obtain the spatial context coding feature sequence.

[0465] Specifically, the dynamic spatial context encoding method includes the following sub-steps:

[0466] S401, Obtain the multidimensional input features of all monitoring stations within the target area at the current discrete time step. The multidimensional input features include at least the air quality observation data and meteorological observation data of the monitoring stations, as well as the location forecast information extracted from the wide-area background field.

[0467] S402, parameterizes the input features of each monitoring station through a first transformation, a second transformation, and a third transformation, generating a source representation, a content representation, and a query representation, respectively. The source representation provides information that can be queried by other stations, the content representation carries the core features to be aggregated, and the query representation expresses the query requirements of the target station.

[0468] S403, based on an attention mechanism, dynamically generates the spatial aggregation weight of the target station at the current time step by interacting with the query representation of each target station and the source representation of all monitoring stations. The spatial aggregation weight reflects the spatial correlation strength between the target station and other monitoring stations and can change in real time with meteorological conditions such as wind direction and speed.

[0469] S404, the content representations of all monitoring stations are weighted and aggregated using the spatial aggregation weights to generate the spatial context code of the target station at the current time step. The above process is repeated for all monitoring stations within the target area to obtain the spatial context codes of all stations at the current time step.

[0470] S405, repeat steps S401 to S404 for multiple consecutive historical time steps, and combine the spatial context codes of all monitoring stations generated in each time step in chronological order to form a three-dimensional tensor with the shape of (number of monitoring stations, length of historical time step, and spatial code dimension), which is the spatial context code feature sequence.

[0471] It should be noted that the specific implementation of the dynamic spatial context encoding method can be found in the detailed description of Embodiment 3 of this specification, and will not be repeated here. Those skilled in the art should understand that the generation method of the spatial context encoding feature sequence is not limited to the method described in Embodiment 3; any feature sequence capable of characterizing the dynamic spatial interaction between monitoring stations can be applied to this embodiment.

[0472] S410, the spatial context encoded feature sequence is input into the pluggable time series prediction module, and the predicted concentration values ​​of atmospheric environmental factors for each monitoring station in the future time step are output.

[0473] The pluggable timing prediction module is configured with a Mixture of Experts (MoE) architecture and features standardized input and output interfaces. The design and implementation of this module are described in detail below.

[0474] (a) Standardized Interface Design

[0475] The pluggable timing prediction module adopts standardized input and output interfaces to ensure the module's independence and replaceability.

[0476] Standardized input interface: Defines a fixed format for the data received by the module. In this embodiment, the input is a four-dimensional tensor of shape (batch_size, num_stations, history_len, feature_dim), where batch_size is the batch size, num_stations is the number of monitoring stations, history_len is the length of the historical time step, and feature_dim is the dimension of the context-aware feature.

[0477] Standardized output interface: Defines a fixed format for the module's output prediction results. In this embodiment, the output is a four-dimensional tensor of shape (batch_size, num_stations, pred_len, target_dim), where pred_len is the prediction time step length and target_dim is the dimension of the prediction target (for example, if predicting only PM2.5, then target_dim=1; if predicting multiple pollutants such as PM2.5, PM10, and O3 simultaneously, then target_dim equals the number of atmospheric environmental factors).

[0478] This standardized interface design enables the pluggable timing prediction module to operate independently of the upstream spatial context coding module and the downstream application module. Any timing prediction model instance conforming to this interface specification can be integrated as an expert into the hybrid expert architecture of this invention, achieving flexible, plug-and-play expansion. For example, when a more advanced timing prediction algorithm emerges, it can be encapsulated as an expert subunit conforming to the interface specification to replace the original expert without modifying other parts of the system.

[0479] (II) Composition of Hybrid Expert Architecture

[0480] The hybrid expert architecture includes the following core components:

[0481] K expert prediction subunits: Each expert prediction subunit is an independent temporal prediction model instance used for parallel processing of the input spatial context-encoded feature sequence. To enhance the model's adaptability to different site types and pollution patterns, the expert prediction subunits can be configured with different types of temporal prediction models, such as Long Short-Term Memory (LSTM), Gated Recurrent Unit (GRU), Transformer, Informer, PatchTST, iTransformer, etc., or any combination of the above models. Each expert subunit is independent of the others and can be specialized for modeling specific patterns.

[0482] Gated network: Used to dynamically generate the weight distribution of K experts based on input-related information. In this embodiment, the input to the gated network includes at least one of the following:

[0483] Site identification information: The unique identifier of each monitoring site is mapped to a learnable low-dimensional dense vector through an embedding layer. This vector represents the static attributes of the site (such as geographical location, site type, etc.).

[0484] Derived features of spatial context-coded feature sequences: Features obtained by statistically aggregating the spatial context-coded feature sequences, including the mean, variance, and trend slope of each monitoring station's historical time steps, are used to characterize the overall dynamic characteristics of the sequence.

[0485] The gated network preferably adopts a self-attention gated network structure, and its working principle is as follows:

[0486] The identifier of each monitoring station is embedded as the query.

[0487] The internal states of the encoders of each expert-predicted sub-unit (such as the hidden state of the last time step) are used as keys and values.

[0488] The attention calculation layer dynamically evaluates the degree of fit of each expert prediction sub-unit to the current input data, and generates attention weights as the expert weight distribution.

[0489] Specifically, the gated network uses a self-attention mechanism to dynamically generate expert weights. Its core principle is to ensure that the weights depend not only on the static attributes of the sites but also on the real-time responses of each expert to the current input data. The specific generation process is as follows:

[0490] (1) Input preparation

[0491] For each monitoring station i, the gating network receives the following inputs:

[0492] Query vector: Identifier embedding q for site i i ∈R E1The site ID is mapped to a dense vector through a learnable embedding layer, representing the static attributes of the site (such as geographic location and site type), and E1 is the dimension of the query vector.

[0493] Key vector set: For each expert prediction subunit r, obtain its internal state h encoded by the current input sequence. ir ∈R E2 For example, the last hidden state of an LSTM or the [CLS] labeled output of a Transformer. This internal state reflects the expert prediction subunit r's understanding of the current input data at site i and the result of feature extraction, where E2 is the dimension of the key vector.

[0494] Value vector set: Expert context representations v that correspond one-to-one with key vectors ir ∈R E3 The value vector can be obtained directly using the key vector or through a linear transformation, where E3 is the dimension of the value vector. In the subsequent weighted fusion step, the value vector is the prediction output of each expert for site i itself, or the result of its linear transformation.

[0495] (2) Calculation of attention score

[0496] For site i, calculate the similarity score between its query vector and each expert key vector:

[0497] , or S ir =MLP([q i h ir ]).

[0498] Among them, S ir Let be the similarity score between the query vector of site i and the r-th expert key vector. is a scaling factor used to stabilize the gradient; MLP is an optional nonlinear mapping. The score S ir This characterizes the degree of matching between the static attributes of site i and the understanding of the current data by the expert prediction subunit r.

[0499] (3) Weight normalization

[0500] Softmax normalization is applied to the scores of all experts to obtain the final expert weight distribution:

[0501] w ir =exp(S ir ) / (exp(S i1 )+……+exp(S iK )).

[0502] The weight distribution satisfies that the sum of all weights equals 1. Furthermore, as the input data changes dynamically, when the features of the input sequence change, the internal state of each expert changes accordingly, leading to a corresponding adjustment in attention scores and weight distribution.

[0503] (4) The dynamic nature of technology

[0504] The dynamism of the above mechanism is reflected on two levels:

[0505] Data-driven: The weights are determined not only by the site ID, but also by the real-time encoding results of the experts on the current time series. For the same site, the sequence features input at different time steps are different, the internal states of the experts are different, and the resulting weight distribution also changes accordingly.

[0506] Interactive perception: The interactive computation of query vector and key vector enables the gating network to explore each expert's understanding of the current input and select the expert who best fits the current context for weighting.

[0507] (5) Comparison with static gating

[0508] Static gating determines weights solely based on site IDs or fixed rules, and these weights remain constant regardless of input data, making them suitable for scenarios with highly stable site patterns. In contrast, the self-attention gating of this invention combines site IDs with expert understanding of the current input to determine weights. These weights can dynamically change with the input sequence, making it suitable for scenarios where site patterns change over time and require real-time adaptation.

[0509] Through the above design, the gating network can achieve true dynamic weight allocation: when traffic stations experience morning rush hour pollution, experts skilled at capturing steep increase patterns are assigned high weights; when background stations exhibit stable fluctuations, experts skilled at modeling periodic patterns are activated. This dynamic adaptive capability is the core innovation of this invention.

[0510] It should be noted that the specific parameters such as the similarity measurement method (dot product or MLP), vector dimension, and activation function in the above calculation process can be adjusted by those skilled in the art according to actual application needs, and all of these fall within the protection scope of this invention.

[0511] This self-attention gating mechanism not only considers the inherent attributes of the site, but also incorporates the expert's understanding of the current input data, enabling highly adaptive and context-aware predictive fusion.

[0512] (III) Weighted fusion output

[0513] The outputs of the K expert prediction sub-units are weighted and summed according to the weight distribution generated by the gating network to obtain the final predicted values ​​of atmospheric environmental factor concentrations. For each monitoring station i, the prediction result is... A weighted combination of the prediction outputs from each expert prediction subunit: .in, The prediction output of the r-th expert prediction subunit for station i is the value vector v. ir This weighted fusion process preserves the specialized contributions of each expert, enabling the model to dynamically select the most suitable prediction mode for different sites.

[0514] It should be noted that the specific parameters in the above calculation process, such as the similarity measurement method (dot product or MLP), vector dimensions E1, E2, E3, and activation function, can be adjusted by those skilled in the art according to actual application needs, and all fall within the protection scope of this invention. This self-attention gating mechanism not only considers the inherent attributes of the site but also combines the expert's understanding of the current input data, enabling highly adaptive and context-aware predictive fusion.

[0515] The hybrid expert time series prediction method in this embodiment has the following significant advantages through standardized interfaces and modular design:

[0516] The algorithm is pluggable: expert sub-units can be flexibly replaced or added or deleted according to actual needs without reconstructing the entire model;

[0517] Site adaptation: The gating network dynamically assigns weights based on site identifiers and input features, enabling different sites to obtain predictions that match their local characteristics;

[0518] Easy to expand: When a new time series prediction algorithm needs to be introduced, it can be seamlessly integrated simply by encapsulating it as an expert sub-unit that conforms to the interface specification.

[0519] Furthermore, the K expert prediction sub-units are collaboratively optimized through differentiated training objectives, enabling each expert to learn complementary time-series prediction patterns, avoiding expert homogenization, and improving overall prediction performance.

[0520] In standard hybrid expert architectures, if all expert prediction subunits use the same training objective and data distribution, it can easily lead to convergence in the patterns learned by the experts, i.e., expert homogenization, thereby weakening the advantages of hybrid expert architectures. To solve this problem, this embodiment introduces differentiated training strategies, specifically including the following two implementation methods, which can be used individually or in combination:

[0521] (i) Differentiated training based on sample allocation

[0522] Different subsets of training samples are assigned to different expert prediction subunits, allowing each expert to focus on specific types of time series patterns. Specifically:

[0523] Assignment by Site Type: Based on the static attributes of the monitoring sites (such as traffic stations, industrial zone stations, background stations, residential zone stations, etc.), the training samples are grouped by site type, and each group of samples is used to train the corresponding expert sub-unit. For example, samples from traffic stations are mainly assigned to LSTM experts who are good at capturing steep rise and fall patterns, while samples from background stations are mainly assigned to Transformer experts who are good at modeling stable periodic patterns.

[0524] Assignment by time pattern: Based on the time characteristics of the training samples (such as morning and evening peak hours, low-emission periods at night, special patterns during holidays, etc.), samples from different time periods are assigned to different experts. For example, samples from peak hours are assigned to Gated Recurrent Unit (GRU) experts who are sensitive to mutations, while samples from stable nighttime periods are assigned to Informer experts who have strong capabilities in modeling long-term dependencies.

[0525] Allocation by pollution event type: Based on the types of historical pollution events (such as sandstorms, secondary formations, straw burning, etc.), samples of different event types are allocated to specialized expert sub-units, enabling each expert to learn the evolution patterns of specific pollution processes.

[0526] During training, each expert updates parameters only on their assigned subset of samples, thus forming a specialized division of labor. It should be noted that sample allocation can be "hard allocation" (each sample is assigned to only one expert) or "soft allocation" (each sample is assigned to multiple experts with a certain probability), and those skilled in the art can choose according to actual needs.

[0527] (ii) Collaborative optimization based on diversity penalty terms

[0528] An expert diversity penalty term is introduced into the loss function to encourage experts to learn complementary prediction patterns. Specifically, the total loss function L... total Defined as:

[0529] L total =L prediction +λ×L diversity .

[0530] Among them, L prediction To predict loss, it is used to measure the difference between the predicted value and the actual value; L diversity λ is a diversity penalty term used to encourage differences among experts; λ is a balance coefficient used to adjust the strength of the diversity penalty, which needs to be adjusted appropriately: too large a value may cause experts to focus only on differences and ignore prediction accuracy, while too small a value will not effectively promote diversity.

[0531] Diversity penalties may take at least one of the following forms:

[0532] (1) Output diversity penalty:

[0533] .

[0534] and These represent the predicted output values ​​of the g1-th and g2-th expert prediction subunits for the same input sample, respectively.

[0535] This penalty encourages different experts to make their predictions for the same input sample as different as possible, so that each expert learns a differentiated mapping function, thereby avoiding the convergence of prediction patterns among experts (i.e., expert homogenization), and enabling them to provide complementary information in subsequent gated weighted fusion.

[0536] (2) Representation-based diversity penalty:

[0537] .

[0538] Among them, h g1 and h g2 These represent the internal states of expert prediction subunits g1 and g2 after encoding the input sequence (e.g., the last hidden state of an LSTM or the [CLS] flag output of a Transformer). sim(h) g1 h g2 The similarity between the two is denoted as cosine similarity. This penalty term encourages different experts to make their internal representations as dissimilar as possible, thereby prompting each expert to understand the input data from different perspectives.

[0539] (3) Gradient-based diversity penalty:

[0540] .

[0541] This represents the gradient direction of the loss function with respect to the parameters of the g1 expert prediction subunits under the current training samples. This represents the gradient direction of the loss function with respect to the parameters of the g2 expert prediction subunits under the current training samples. This represents the calculation of the squared Euclidean distance between two gradient vectors, which measures the degree of difference in the update directions of the two expert parameters.

[0542] This penalty term encourages different experts to adopt different gradient directions by maximizing the distance (i.e. minimizing its negative value), thereby guiding each expert to update parameters in different optimization directions and learn complementary prediction patterns.

[0543] (III) Combined Training Strategies

[0544] Sample allocation can be either "hard allocation" (each sample is assigned to only one expert, suitable for scenarios with high differentiation between site types / contamination patterns) or "soft allocation" (each sample is assigned to multiple experts with a certain probability, suitable for complex scenarios with cross-fusion of contamination patterns), and those skilled in the art can choose according to the characteristics of the dataset.

[0545] The two methods described above can be combined to form a dual differentiation mechanism of sample allocation and diversity penalty: first, hard allocation allows experts to initially form a division of labor, and then the diversity penalty term further strengthens the complementarity among experts. For example, in the early stages of training, a sample allocation strategy is used to enable experts to quickly converge to different pattern spaces; in the later stages of training, a diversity penalty term is introduced to finely adjust the boundaries between experts and avoid pattern overlap caused by blurred boundaries.

[0546] Through the above-described differentiated training strategy, this embodiment achieves the following beneficial effects:

[0547] Enhanced expert complementarity: Each expert learns differentiated time series forecasting patterns, avoiding homogenization problems and fully leveraging the integrated advantages of the hybrid expert architecture;

[0548] Overall prediction accuracy improved: When faced with complex and ever-changing pollution processes, the gating network can dynamically select the most suitable expert combination based on the input data, thereby improving the overall prediction accuracy.

[0549] Enhanced robustness: Due to the complementarity among experts, even if an individual expert makes a large prediction error in a specific situation, other experts can still compensate through weighted fusion, thereby enhancing the robustness of the model.

[0550] It should be noted that the above-mentioned sample allocation method, the specific form of the diversity penalty term, and the value of the balance coefficient λ can be adjusted by those skilled in the art according to the actual application scenario and the characteristics of the dataset, and all of these fall within the protection scope of this invention.

[0551] Furthermore, the method also includes an end-to-end training step: by minimizing the loss function between the predicted value and the true value, the parameters of the K expert prediction subunits and the gating network are synchronously updated using the backpropagation algorithm.

[0552] To achieve collaborative optimization of the hybrid expert architecture, this embodiment adopts an end-to-end training approach, treating the entire pluggable temporal prediction module (including K expert prediction sub-units and a gating network) as a single model, and synchronously updating all learnable parameters through a unified loss function and backpropagation algorithm. Specifically:

[0553] (I) Loss Function Design

[0554] The total loss function consists of the prediction loss and an optional diversity penalty term:

[0555] L total =L prediction +L reg .

[0556] Among them, L prediction To predict loss, it measures the difference between the model's predicted values ​​and the actual observed values. Depending on the nature of the prediction task, mean squared error (MSE), mean absolute error (MAE), or negative log-likelihood (NLL) can be selected. For deterministic predictions, MSE is preferred; for probabilistic predictions, NLL is preferred. reg Regularization terms are optional components that can be introduced with different constraints depending on specific needs. For example, to encourage experts to learn complementary prediction patterns, the diversity penalty terms mentioned above (such as output-based Euclidean distance penalty, representation-based cosine similarity penalty, or gradient-based diversity penalty) can be used; to prevent overfitting, L1 or L2 regularization can be used; and to improve prediction smoothness, temporal consistency constraints can be introduced. Those skilled in the art can select appropriate regularization terms and their combinations according to actual needs.

[0557] (II) Forward Propagation and Loss Calculation

[0558] During training, for each training sample (i.e., a sequence of historical spatial context encoded features from a batch of monitoring stations):

[0559] The input is simultaneously fed to K expert prediction sub-units, and each expert performs forward calculations independently to obtain their own prediction outputs.

[0560] The gating network calculates the weight distribution for each expert based on site identifiers and / or derived features of the input sequence;

[0561] The predictions from each expert are weighted and fused to obtain the final prediction value.

[0562] The prediction loss is calculated based on the final predicted value and the actual value, and is also calculated based on the selected regularization term;

[0563] Calculate the total loss.

[0564] (III) Backpropagation and parameter update

[0565] The gradient of the total loss function with respect to the parameters of each expert prediction subunit and the gating network parameters are calculated using the backpropagation algorithm. The gradient reflects the direction and extent of each parameter's influence on the loss function and serves as the basis for parameter adjustment.

[0566] Subsequently, gradient descent optimizers (such as Adam, SGD, etc.) are used to simultaneously update all learnable parameters. The update rule follows the standard gradient descent principle: each parameter is adjusted along its negative gradient direction with a step size controlled by the learning rate. Specifically, for the parameters of each expert prediction subunit and the parameters of the gated network, they are iteratively updated according to the formula "new parameter = old parameter - learning rate × gradient".

[0567] The entire training process involves repeated forward propagation, loss calculation, backpropagation, and parameter updates until the model converges or reaches a preset stopping condition. Through this end-to-end joint optimization, the expert prediction subunits and the gating network can learn collaboratively to improve prediction performance. Through end-to-end joint training, the expert prediction subunits and the gating network can optimize together: the gating network learns how to dynamically allocate expert weights based on the input, while each expert gradually develops a specialized division of labor through interaction with the gating network. This joint optimization approach enables the entire hybrid expert architecture to adaptively capture pollution time-series patterns at different sites and time periods, thereby improving overall prediction performance.

[0568] It should be noted that the specific form of the regularization term in the above loss function, the value of the balancing coefficient (if needed), the selection of the optimizer, etc., can be adjusted by those skilled in the art according to the actual application scenario, and all of these adjustments fall within the protection scope of this invention. If no regularization term (i.e., L) is enabled... reg If the error is 0, the training process degenerates into the standard end-to-end prediction loss minimization, which still falls within the protection scope of this invention.

[0569] In this embodiment, the predicted atmospheric environmental factor concentration is used for probabilistic downscaling and uncertainty estimation of atmospheric environmental factor concentration at the station level. It is input into the parameterized network as part of the multimodal input information, and the parameterized network outputs the probability distribution parameters of any predicted location.

[0570] Furthermore, the method provided in this embodiment also includes a probabilistic downscaling step:

[0571] The predicted concentrations of atmospheric environmental factors at each monitoring station in future time steps are used as part of the multimodal input information and are input into the parameterized network along with the wide-area background field data and static geographic attribute data.

[0572] Given the geographic coordinates of any predicted location within the target area, the parameterized network outputs a set of parameters representing the probability distribution of atmospheric environmental factor concentrations at the predicted location in future time steps.

[0573] Based on the probability distribution parameter set, a complete probability density function for the concentration of the atmospheric environmental factors is constructed.

[0574] The specific implementation of this step will be described in detail in Example 6.

[0575] Based on the same inventive concept, this embodiment also provides a hybrid expert time series prediction device for local site prediction, comprising:

[0576] The input module is used to obtain spatial context-coded feature sequences that characterize the spatial interactions of each monitoring station.

[0577] The prediction module includes a pluggable time-series prediction module, which has a standardized input interface and a standardized output interface, and is configured to receive the spatial context-coded feature sequence and output the predicted atmospheric environmental factor concentration values ​​for each monitoring station at future time steps.

[0578] The pluggable timing prediction module is configured as a hybrid expert architecture, which includes:

[0579] There are K expert prediction subunits, each of which is an independent temporal prediction model instance, used to process the spatial context encoded feature sequence in parallel.

[0580] A gating network is used to dynamically generate the weight distribution of K expert prediction sub-units based on site identification information and / or derived features of the spatial context encoded feature sequence.

[0581] The fusion unit is used to weight and fuse the outputs of the K expert prediction subunits according to the weight distribution to obtain the predicted values ​​of the atmospheric environmental factor concentrations.

[0582] The local site spatiotemporal prediction model (local expert) of the present invention is composed of the dynamic spatial context encoding module described in Embodiment 3 and the pluggable temporal prediction module described in Embodiment 4. Its overall workflow and collaboration mechanism are as follows:

[0583] I. Overall Workflow of Local Model

[0584] 1. Data Input Layer

[0585] The local model receives three types of input data:

[0586] Site observation data: historical atmospheric environmental factor concentration sequences (such as PM2.5, PM10, O3, etc.) and local meteorological observation data (such as temperature, humidity, wind speed, wind direction) of each monitoring station.

[0587] Wide-area background field data: The grid values ​​and neighborhood information corresponding to each station location extracted from the forecast field output in Example 2 or Example 5 are used as the external driving background;

[0588] Auxiliary information includes station identification information (used to distinguish the static attributes of different stations), static geographic attribute data (such as elevation, land use type), time characteristics (such as hour, weekday, and holiday identifiers), and optional future weather forecast data.

[0589] After preprocessing, the above data forms an input tensor of shape (num_stations, history_len, input_dim), where history_len is the length of the historical time step.

[0590] 2. Spatial coding layer (corresponding to Example 3)

[0591] The input data is first fed into the dynamic spatial context encoding module. This module performs the following operations independently at each discrete time step:

[0592] For each monitoring station, the input features are parametrically transformed by a first, second, and third transformation to generate source representation, content representation, and query representation, respectively.

[0593] Based on the attention mechanism, the query representation of the target site and the source representation of all sites are used to perform interactive calculations to dynamically generate spatial aggregation weights. These weights can change in real time with meteorological conditions such as wind direction and wind speed, reflecting the time-varying pollution transmission relationship between sites.

[0594] The spatial aggregation weights are used to perform weighted aggregation of the content representations of all sites to generate the spatial context encoding of the target site at the current time step.

[0595] The above process is repeated for multiple consecutive historical time steps, and the spatial context codes generated at each time step are combined in chronological order to obtain a spatially encoded feature sequence of shape (num_stations, history_len, spatial_dim). This sequence fully preserves the dynamic spatial interaction information of each station at each historical time step.

[0596] 3. Time Series Prediction Layer (corresponding to Example 4)

[0597] The spatially encoded feature sequence is input to the pluggable temporal prediction module. This module has standardized input and output interfaces and supports two operating modes:

[0598] Single model mode: The input sequence is directly processed by the selected base time series model (such as LSTM, Transformer, etc.) to output the predicted concentration values ​​of atmospheric environmental factors for each monitoring station in the future time step.

[0599] Hybrid Expert (MoE) Mode: Data is simultaneously fed into K expert prediction sub-units and a gating network. Each expert prediction sub-unit is an independent time-series prediction model instance, processing the input sequence in parallel. The gating network dynamically generates the weight distribution of the K experts based on site identification information and / or derived features of the input sequence. The prediction outputs of each expert are weighted and fused according to the weight distribution to obtain the final prediction result. This mode can perform specialized modeling for different site types or different pollution processes, achieving adaptive prediction.

[0600] 4. Output Layer

[0601] The time-series prediction layer outputs a prediction tensor of shape (num_stations, pred_len, target_dim), where pred_len is the prediction time step length and target_dim is the dimension of the prediction target (e.g., the number of atmospheric environmental factors). For each monitoring station, a refined air quality prediction value for a specified future time step can be obtained.

[0602] 5. End-to-end training mechanism

[0603] The entire local model is jointly trained end-to-end. A suitable loss function (e.g., mean squared error (MSE), mean absolute error (MAE), or negative log-likelihood (NLL)) is defined to measure the difference between the predicted and true values. The gradient of the loss function with respect to all learnable parameters in the model is calculated using the backpropagation algorithm, including:

[0604] The parameters of the first, second, and third transforms in the spatial coding module;

[0605] The parameters of each expert prediction subunit in the time series prediction module (the parameters of the basic time series model in the single model mode).

[0606] Parameters of the gating network (in MoE mode).

[0607] Gradient descent optimizers (such as Adam and SGD) are used to update all parameters synchronously and iteratively optimize until the model converges. This end-to-end joint training approach enables the spatial encoding module and the temporal prediction module to optimize collaboratively: the spatial encoding module learns to generate spatial features that are beneficial to temporal prediction, while the temporal prediction module improves prediction accuracy with the support of spatial features, thereby achieving better overall prediction performance.

[0608] II. Characteristics of Local Model Technology

[0609] 1. Expression of attention weights influencing site relationships

[0610] Through the attention mechanism in the dynamic spatial context encoding module, the attention weights generated by the model at each time step reflect the spatial influence relationships between sites. These weights change with meteorological conditions such as wind direction and speed, and can be used for visualization analysis to help environmental managers identify the sources, transport paths, and diffusion trends of pollutants, providing a reference for pollution source tracing and emergency decision-making.

[0611] 2. Co-modeling of spatiotemporal coupling characteristics

[0612] The local model captures two types of spatiotemporal features—"dynamic spatial interaction between sites" and "temporal evolution patterns within sites"—through the collaborative design of two modules: spatial encoding and temporal prediction. The spatial encoding module incorporates the wide-area transmission background into the local feature representation, while the temporal prediction module employs a hybrid expert architecture to adapt to local contamination patterns at different sites. These mechanisms work together to predict complex contamination processes.

[0613] 3. Module decoupling and the generalization capability of hybrid expert architecture

[0614] Modular design decouples spatial coding from temporal prediction, allowing each module to be optimized or updated independently. The hybrid expert architecture dynamically selects or combines multiple expert models through a gating network, enabling it to adapt to site types or data patterns not found in the training set and exhibiting low overfitting to specific training sets. Furthermore, the prediction results from multiple experts mutually correct each other during the weighted fusion process, maintaining a certain level of performance even under abnormal conditions.

[0615] 4. The indicative role of gating weights in pollution driving factors

[0616] The gating weight allocation results in the hybrid expert architecture can be used to analyze the pollution drivers at different sites. By examining the types of experts activated at a specific site and their weight changes, the dominant factors of the site's pollution patterns (such as traffic emissions, industrial transmission, secondary generation, etc.) can be inferred, providing directional references for refined pollution prevention and control.

[0617] 5. Standardized interfaces and component replaceability

[0618] The time-series prediction module adopts a standardized interface design and is pluggable. In practical applications, the internal prediction algorithm can be replaced or upgraded as needed without modifying the spatial coding module or other system components. This design facilitates technology iteration and supports rapid adaptation to different regions, pollutants, or application scenarios.

[0619] Example 5

[0620] This embodiment relates to a probabilistic forecasting method for wide-area atmospheric environmental factors, which belongs to the uncertainty quantification link of the wide-area model in the dual-expert fusion prediction architecture of this invention. It is used to output the covariance matrix based on deterministic forecasting, and provide uncertainty priors for subsequent local probabilistic prediction (Example 6).

[0621] In wide-area grid forecasting, uncertainties arise from extremely complex sources, primarily including inherent errors in the input weather forecast field, approximation errors in the model itself, and the inherent chaotic characteristics of the atmospheric system. However, existing wide-area forecasting models suffer from the following technical shortcomings:

[0622] The contradiction between computational cost and timeliness: Traditional numerical weather prediction models (such as WRF-Chem) are based on solving partial differential equations of physical mechanisms, which has huge computational overhead. A single high-resolution forecast often takes several hours, making it difficult to meet the operational needs of real-time updates.

[0623] Insufficient ability to model spatial long-range dependence and uncertainty: Existing machine learning models (such as convolutional neural networks) are limited by local receptive fields and have difficulty effectively capturing the long-distance transport of pollutants at the kilometer level; at the same time, mainstream deep learning methods (such as Google GraphCast and Huawei Pangu meteorological big data model) mainly focus on deterministic forecasts, outputting a single atmospheric state field at a future moment, which cannot quantify the uncertainty of forecast results and makes it difficult to assess the reliability of forecast results for downstream tasks (such as site downscaling and extreme event warning).

[0624] Challenges in fusing multi-source heterogeneous data: Ground observation fields and vertical profile data differ significantly in spatial resolution, sampling methods, and physical meaning. Existing technologies struggle to achieve effective joint modeling of the two, resulting in insufficient utilization of atmospheric structure information in the vertical direction.

[0625] Lack of physical consistency: Pure data-driven models lack explicit constraints on atmospheric physical laws, and the forecast results may violate basic physical laws such as the conservation of mass and momentum, reducing the reliability of the forecast results.

[0626] To address the aforementioned technical problems, this embodiment provides a probabilistic forecasting method for the concentration of atmospheric environmental factors over a wide area. The core concept of this method is as follows:

[0627] Efficient spatiotemporal modeling: Employing a three-dimensional spatiotemporal attention mechanism, it simultaneously captures the dynamic correlation of the atmosphere in horizontal space, vertical stratification, and time dimension, achieving accurate modeling of pollution transmission at the kilometer level;

[0628] Uncertainty quantification: The covariance matrix of the forecast uncertainty is output by the covariance prediction head, and together with the deterministic forecast field output by the mean prediction head, it defines the multivariate Gaussian distribution of the atmospheric state at future time, providing a pixel-by-pixel probability confidence index.

[0629] Multi-source data fusion: Design a vertical field dynamic projection module to map non-uniformly distributed vertical observation data to a unified grid, realizing joint modeling of ground observation field and vertical profile data;

[0630] Physical consistency guarantee: Physical constraints such as mass conservation are introduced into the loss function, so that the deep learning model learns to follow the physical laws of atmospheric pollutant transport during the training process, thereby enhancing the physical rationality of the forecast results.

[0631] This embodiment provides a method for probabilistic forecasting and covariance modeling of wide-area atmospheric environmental factors. The method is implemented through a pre-constructed neural network model. The neural network model includes a three-dimensional spatiotemporal attention encoder, a mean prediction head, and a covariance prediction head, which learns the mapping relationship from historical observations to future probability distributions through end-to-end training. The following is in conjunction with the appendix... Figure 5 The specific implementation steps of this method are explained in detail.

[0632] like Figure 5 As shown, the method includes the following steps:

[0633] S500, obtains the input spacetime cube of historical time steps.

[0634] This step is the data preparation and preprocessing stage for wide-area probabilistic forecasting. Its core task is to collect and organize multi-source heterogeneous atmospheric observation data, construct an input spatiotemporal cube that can fully describe the atmospheric state of the target area, and provide standardized input for subsequent three-dimensional spatiotemporal attention coding.

[0635] (a) The composition of the input spacetime cube

[0636] The input spacetime cube is a four-dimensional tensor structure encompassing four dimensions: time, spatial height, spatial width, and atmospheric parameter channels. Specifically, it includes the following information:

[0637] Time dimension: Includes data from multiple consecutive historical time steps. In this embodiment, it is preferable to use data from the current time T and the previous time T-1 to form an input sequence with a time span of 2, in order to capture the short-term temporal evolution characteristics of atmospheric states. It should be noted that the number of time steps can be adjusted according to the actual prediction task requirements, for example, using more historical time steps to enhance the model's ability to perceive temporal patterns.

[0638] Spatial Dimension: A regular grid covering the target prediction area, defined by spatial height H and spatial width W. The grid resolution can be set according to the application scenario, such as a 0.25°×0.25° latitude and longitude grid or a 3km×3km projected grid. Each grid cell represents a spatial location, and the atmospheric conditions within that grid will be predicted subsequently.

[0639] Parameter channel dimensions: Includes various atmospheric parameters, divided into two parts: ground observation field and vertical observation field.

[0640] Ground observation field parameters mainly include near-surface atmospheric environmental factor concentration data (such as PM2.5, PM10, SO2, NO2, O3, CO, etc.) and meteorological element data (such as temperature, humidity, air pressure, wind speed, wind direction, precipitation, etc.). These parameters are usually derived from reanalysis data (such as ERA5, GFS / FNL) or gridded interpolation products from ground monitoring stations.

[0641] Vertical observation field parameters include atmospheric composition or meteorological data at different altitudes, such as ozone concentration, temperature, humidity, and wind speed at multiple isobaric surfaces (e.g., 1000 hPa, 850 hPa, 700 hPa, 500 hPa). Vertical observation data plays a crucial role in understanding the vertical transport of pollutants, boundary layer evolution, and regional transport.

[0642] (II) Data Sources and Preprocessing

[0643] The data involved in the input spatiotemporal cube comes from multiple sources and needs to undergo a unified preprocessing procedure to eliminate the differences between the multiple data sources:

[0644] Data sources primarily include global or regional numerical model reanalysis data (such as ERA5 and CAMS reanalysis), operational meteorological forecast products (such as GFS and ECMWF forecast fields), satellite remote sensing inversion products, and observational data from ground monitoring stations. For real-time forecasting tasks, the corresponding forecast field or analysis field can be used as input.

[0645] Spatiotemporal alignment: Since different data sources may have different spatial resolutions, temporal resolutions, and projection methods, it is necessary to unify all data to the same spatial grid and temporal frequency through interpolation or resampling methods. For example, meteorological fields from different sources can be mapped to the target grid through bilinear interpolation, and hourly and six-hourly data can be processed through temporal interpolation or alignment.

[0646] Missing value handling: For observational data with missing values, preprocessing methods such as nearest-neighbor spatiotemporal interpolation, climatological substitution, or model prediction imputation are used to ensure the integrity of the input data.

[0647] (III) Processing of Vertical Observation Data: Steps for Dynamic Projection of Vertical Field

[0648] Since vertical observation data is typically distributed across non-uniform altitude or pressure layers, and the definitions of vertical layers from different sources may differ, directly stacking them as channels would make it difficult for the model to effectively utilize vertical structure information. To address this issue, this embodiment introduces a dynamic vertical field projection step, mapping the non-uniformly distributed vertical observation data onto a unified vertical grid.

[0649] 1. Definition of raw vertical observation data

[0650] For each spatial grid point, the acquired raw vertical observation data can be represented as:

[0651] The number of original vertical layers L orig It may change with spatial location;

[0652] The layer identifier p corresponding to each original layer L (such as air pressure, altitude);

[0653] The observed feature vector f on each original layer L ∈R Cvert Where Cvert is the number of vertical observation parameter channels.

[0654] 2. Setting a uniform vertical grid

[0655] Pre-define a set of uniform vertical layer target positions {L1, ..., L...} u , ..., L M}, where L u Let u be the target location for the u-th vertical layer, where u ranges from 1 to M, and M is the number of uniform vertical layers (e.g., 15 isobaric surfaces). These target locations are spatially consistent across all grid points, allowing subsequent processing to be performed within a unified vertical reference frame.

[0656] 3. Dynamic projection calculation

[0657] For each uniform vertical layer target location L u Its projection features are obtained by weighted aggregation of features from all original vertical layers. The weights are dynamically calculated based on the degree of matching between the original layer positions and the target layer positions. Weight calculation can be performed using one of the following methods:

[0658] Distance-weighted: Weights and |p L -L u | Inversely proportional, the closer the original layer is, the greater its contribution;

[0659] Kernel-weighted: The weights are calculated using a Gaussian kernel function, with the kernel width being a learnable parameter;

[0660] Attention mechanism: The target layer position is encoded as a query vector, the original layer position is encoded as a key vector, and the original layer features are used as value vectors. Weights are calculated through attention.

[0661] The learnable parameters in the aforementioned weight calculation (such as kernel width, positional encoding matrix, and attention query / key transformation matrix) are all jointly optimized with the entire prediction method through end-to-end training. During training, the projection method is automatically adjusted according to the loss function of the downstream prediction task, so that the projected vertical features are most conducive to the final atmospheric environmental factor concentration prediction.

[0662] 4. Projection Result Output

[0663] After performing the above projection calculation on each uniform vertical layer target location, a projected vertical feature with shape (M, Cvert) is obtained. This projected vertical feature is then concatenated with the ground observation field parameters along the channel dimension to form a complete input feature containing vertical information, which is then used by the subsequent 3D spatiotemporal attention encoder.

[0664] The vertical field dynamic projection step is implemented by the vertical field dynamic projection module, which is built into the neural network model. This module is used to map the non-uniformly distributed original vertical observation data onto a preset uniform vertical grid, generating vertical features consistent with the spatial resolution of the ground observation field. (IV) Final form of the input spatiotemporal cube

[0665] After the above data acquisition, preprocessing, and vertical projection, the final spatiotemporal cube used as model input is obtained, with the following dimensions:

[0666] Input tensor shape = (B, Tin, H, W, Ctotal).

[0667] Where: B is the batch size; Tin is the input time step (2 in this embodiment); H and W are the spatial grid height and width; Ctotal=Cbase+Cvert is the total number of parameter channels, Cbase is the number of ground observation field channels, and Cvert is the number of vertical observation field channels after projection.

[0668] The input spatiotemporal cube fully represents the atmospheric state of the target area in the past, including the concentration of atmospheric environmental factors and meteorological elements at the ground layer, and also incorporates vertical distribution information, providing a rich data foundation for subsequent three-dimensional spatiotemporal attention encoding.

[0669] S510, the input spatiotemporal cube is input to an encoder based on a three-dimensional spatiotemporal attention mechanism to extract spatiotemporal feature maps.

[0670] The three-dimensional spatiotemporal attention mechanism used in this step is the same as the three-dimensional window attention mechanism described in S210 of Example 2, that is, by dividing a three-dimensional window, calculating self-attention within the window, and using a window sliding mechanism, it simultaneously captures the dynamic correlation of the atmosphere in the horizontal space, vertical layering, and time dimension. Specifically:

[0671] The input feature map is divided into non-overlapping three-dimensional windows, and the size of each window is the time window size × spatial height window size × spatial width window size;

[0672] Self-attention calculation is performed within each 3D window to perceive the spatiotemporal positional relationship of feature points through relative positional offset.

[0673] Establish cross-window connections through a window sliding mechanism to achieve spatial dissemination of information;

[0674] By employing a hierarchical downsampling structure, the above process is repeated on feature maps at different scales to achieve the extraction of multi-scale spatiotemporal features.

[0675] For details on the specific implementation of the above mechanism, please refer to the relevant description of S210 in Embodiment 2, which will not be repeated here.

[0676] Through processing by a 3D spatiotemporal attention encoder, this step outputs a feature map rich in spatiotemporal information, providing a foundation for subsequent mean and covariance prediction.

[0677] S520, Based on the spatiotemporal feature map, a deterministic forecast mean field for future times is generated through the mean prediction head.

[0678] The goal of this step is to decode the mean concentration field of atmospheric environmental factors for future times from the spatiotemporal feature map extracted from S510, i.e., the deterministic forecast part. The mean prediction head adopts a convolutional neural network structure, taking the spatiotemporal feature map as input and outputting a multidimensional mean field with the same spatial resolution as the input.

[0679] Specifically, the mean prediction head consists of several convolutional layers that progressively map the high-dimensional features output by the encoder to the target number of output channels. First, one or more 3×3 convolutional layers further process the spatiotemporal feature map, fusing local neighborhood information. Then, 1×1 convolutional layers adjust the number of channels to Cbase, which is the number of atmospheric environmental factors (such as PM2.5, PM10, O3, etc.) in the ground observation field. Batch normalization and nonlinear activation functions (such as ReLU, GELU) can be inserted between convolutional layers to enhance expressive power.

[0680] The deterministic forecast mean field μc∈R output by the mean prediction head H×W×CbaseIt is a multidimensional tensor with the same size as the input spatial grid, where each spatial location corresponds to a Cbase-dimensional vector representing the predicted concentration of each pollutant at that location at future times. This deterministic forecast mean field characterizes the most likely evolution of the atmospheric state and serves as the central estimate of the subsequent probability distribution.

[0681] S530, Based on the spatiotemporal feature map, a covariance matrix for forecasting uncertainty is generated through the covariance prediction head.

[0682] The goal of this step is to quantify the uncertainty of the forecast results and output a covariance matrix describing the joint probability distribution of atmospheric states at future times. Parallel to the mean prediction head, the covariance prediction head also takes the spatiotemporal feature map extracted by S510 as input, but outputs a parameterized representation of the covariance matrix.

[0683] (1) Mathematical meaning of covariance matrix

[0684] The covariance matrix Σ∈R (HW)×(HW) This is a high-dimensional matrix where the diagonal elements represent the prediction variance (i.e., the magnitude of uncertainty) of the concentrations of various atmospheric environmental factors at each spatial location, and the off-diagonal elements represent the correlation of prediction errors between different spatial locations or between different pollutants. This matrix, together with the mean field μc, defines the multivariate Gaussian distribution of the atmospheric state at future times:

[0685] P(X T+1 |[X T-1 X T ])=N(μc,Σ).

[0686] (2) The challenge of efficient modeling of covariance matrix

[0687] Directly outputting the complete (HW)×(HW) covariance matrix is ​​computationally and storage-infeasible—for a grid with H=W=100, the matrix dimension is as high as 10. 4 ×10 4 The number of parameters exceeds one hundred million. To address this issue, this embodiment employs a constraint modeling approach that combines block diagonalization with low-rank global correlation, significantly reducing computational complexity while maintaining expressive power.

[0688] (3) Diagonal part of the block: Local correlation modeling

[0689] Diagonal part of the block (diagonal matrix of the block) Σ local This is used to model the correlations between different atmospheric parameters within each spatial location. Specifically, for each spatial grid point, a Cbase×Cbase covariance matrix is ​​independently modeled to represent the correlations between the concentrations of different atmospheric environmental factors at that location. The local covariance matrices of all spatial locations form a block diagonal matrix:

[0690] Σ local =diag(Σ1,Σ2,…,Σ v1 , …, Σ H×W ).

[0691] Among them, each Σ v1 Let v1 be the local covariance matrix of the v1-th spatial location, where v1 ranges from 1 to H×W. It can be parameterized to Σ using Cholesky decomposition. v1 =R v1 ×R v1 T , where R v1 It is a lower triangular matrix whose elements are output by the covariance prediction head.

[0692] (4) Low-rank global correlation component: cross-spatial correlation modeling

[0693] Low-rank global correlation component (low-rank global correlation matrix) Σ global This is used to model long-range correlations between different spatial locations, such as the spatial association of the same pollutant between different grid points, or the cross-spatial correlation between different pollutants. To achieve efficient modeling, Σ... global Decomposed into a product of low-rank matrices:

[0694] Σ global =UU T .

[0695] Where, U∈R (H×W)×a As a learnable low-rank matrix, the rank *a* is much smaller than the total number of spatial locations (in this embodiment, *a* = 5). This decomposition method reduces the number of parameters in the covariance matrix from O((H×W)). 2 The cost is reduced to O(H×W×a), making cross-spatial correlation modeling feasible in engineering for the first time.

[0696] (5) Additional diagonal perturbation term

[0697] To ensure the positive definiteness of the entire covariance matrix, a diagonal perturbation term is added to the block diagonal and low-rank global correlation:

[0698] Σ=Σ local ⊕Σ global +b·I.

[0699] Where b is a very small positive number (e.g., 10). -5 I is the identity matrix. This perturbation term ensures that all eigenvalues ​​of the covariance matrix are greater than b, avoiding singularity issues in numerical computation. ⊕ indicates element-wise addition.

[0700] (6) Output of the covariance prediction head

[0701] The covariance prediction header needs to output the parameters required for the above parameterized representation, specifically including:

[0702] For each spatial location, output the parameters of the local covariance matrix (such as the Cholesky factor R). v (the lower triangular part);

[0703] Output all elements of the global low-rank matrix U (a total of H×W×a values).

[0704] The covariance prediction head can adopt a convolutional structure similar to the mean prediction head, and finally use a specific activation function to ensure that the output parameters meet the constraints (such as using softplus or exponential functions to ensure that the variance is positive).

[0705] (7) Coordination with the mean field

[0706] The mean field and covariance matrix together define the multivariate Gaussian distribution of the atmospheric state at future times. During training, the two are jointly optimized using a negative log-likelihood loss function; during inference, based on this distribution, the prediction interval and exceedance probability of any atmospheric environmental factor concentration at any spatial location can be calculated, providing comprehensive uncertainty information for environmental decision-making.

[0707] Furthermore, the method also includes a physical constraint enhancement optimization step, which incorporates atmospheric physical laws as prior knowledge during the parameter optimization process of the neural network model to improve the physical rationality and credibility of the forecast results.

[0708] In the parameter optimization process of the neural network model, the total loss function is composed of data-driven loss and physical constraint loss. The data-driven loss uses the negative log-likelihood function to measure the degree of fit between the predicted multivariate Gaussian distribution and the actual observation; the physical constraint loss is based on the law of conservation of atmospheric mass and is achieved by imposing physical consistency constraints on the predicted mean field.

[0709] Specifically, for each training sample, the pollutant advection flux F = μc × v is calculated using the wind field data at the input time (including horizontal wind speed components u and v) and the mean field μc predicted for future times using the method described above. According to the principle of mass conservation, in a source-sink-free region, the local variation in atmospheric environmental factor concentrations should be balanced with advection transport, meaning the divergence of the advection flux should be as small as possible. Therefore, physical constraints are constructed:

[0710] L physics =||▽·(μc×v)|| 2 .

[0711] The divergence ▽·(μc×v) can be calculated using spatial differencing or automatic differentiation. This constraint term is added to the total loss function in a differentiable form.

[0712] L total =L NLL +λ1×L physics .

[0713] Wherein, λ1 is the equilibrium hyperparameter, used to adjust the strength of the physical constraints.

[0714] During backpropagation, the gradient generated by the physical constraint loss guides the learnable parameters of the method to update in a direction that satisfies mass conservation, making the prediction results more statistically consistent with the physical laws of atmospheric transport. This constraint term can be applied to the entire parameter optimization process, or it can be used as an auxiliary task in the fine-tuning stage after pre-training.

[0715] By introducing physical constraints to enhance optimization, this method effectively suppresses anomalous forecasts that violate basic physical laws (such as abnormal accumulation of atmospheric environmental factor concentrations in the source-free region or contradictions between transport direction and wind field) that may occur in pure data-driven methods while maintaining the high efficiency of data-driven methods. This significantly improves the physical reliability and generalization ability of wide-area probabilistic forecasts.

[0716] Furthermore, the method employs a progressive multi-task training framework to optimize the neural network model in stages, thereby improving the predictive performance of the mean field and covariance matrix, and enhancing the model's ability to perceive vertical structures.

[0717] In the parameter optimization process of the neural network model, a progressive multi-task training strategy is introduced, which divides the training process into three stages, each focusing on a different learning objective. By gradually increasing the task complexity, the instability in the early stage of multi-task joint training is avoided, enabling the model to converge more effectively.

[0718] Phase 1: Mean Pre-training

[0719] This phase only trains the mean prediction head to optimize the accuracy of deterministic predictions. The loss function L... stage1 The mean square error is used to measure the difference between the predicted mean field and the actual observed mean field, μc. true Differences between them:

[0720] L stage1 =MSE(μc,μc true ).

[0721] In this stage, the parameters of the 3D spatiotemporal attention encoder and the mean prediction head are jointly optimized, while the parameters of the covariance prediction head are frozen or temporarily excluded from training. Through this stage of training, the model first learns to capture the main evolutionary patterns of atmospheric states, laying the foundation for subsequent uncertainty modeling.

[0722] Phase Two: Covariance Prediction Training

[0723] In this phase, the parameters of the mean prediction head are frozen, and only the covariance prediction head is trained, focusing on uncertainty modeling. The loss function L... stage2 Using negative log-likelihood, the covariance matrix is ​​calculated based on the mean field of the first-stage pre-training and the output of the current covariance prediction head:

[0724] L stage2 =L NLL (μ fixed ,Σ).

[0725] Where, μ fixed This is the fixed mean field output after the first stage of training. The parameters of the 3D spatiotemporal attention encoder continue to be optimized in this stage to meet the needs of covariance prediction. By freezing the mean prediction head, perturbations to the converged mean field are avoided in the early stages of covariance training, allowing the model to focus on learning the distribution characteristics of the prediction error.

[0726] Phase 3: Joint Fine-tuning and Integration of Auxiliary Tasks

[0727] This phase unfreezes the parameters of all modules and performs end-to-end joint fine-tuning. Total loss function L stage3 It is composed of a weighted average of negative log-likelihood loss and auxiliary task loss:

[0728] L stage3 =L NLL +λ3×L aux .

[0729] Where λ3 is a hyperparameter balancing the weights of the auxiliary tasks. Auxiliary tasks can be added as needed, such as a vertical field prediction task: utilizing features from the encoder's intermediate layers, an additional prediction head outputs meteorological or compositional parameters for the vertical layer, and the mean square error is calculated between this parameter and the actual vertical observations as L. aux This auxiliary task forces the encoder to learn more physically meaningful vertical structural features, thereby improving the performance of the main task (ground field probability forecasting).

[0730] In the third stage of training, a progressive weight scheduling strategy can be adopted. For example, the auxiliary task can be started with a smaller λ3, and then λ3 can be gradually increased after the main task stabilizes to balance the learning of the primary and secondary tasks. Through this progressive multi-task training framework, the model can master deterministic forecasting, uncertainty quantification, and vertical structure perception capabilities step by step, ultimately achieving high-quality, physically reliable wide-area probabilistic forecasts.

[0731] To further improve the deployment efficiency of this method in actual business systems, the following optimization techniques can be adopted:

[0732] 1. Gradient checkpointing technique

[0733] When training large-scale 3D Swing Transformer models, GPU memory usage is a major bottleneck. To address this, gradient checkpointing can be enabled during the forward propagation of 3D Swing Transformer blocks: intermediate activation values ​​are not retained during forward propagation but are recalculated during backward propagation, trading time for space and significantly reducing GPU memory usage. This technique can reduce GPU memory requirements by approximately 50% when training high-resolution (e.g., 1024×1024 grid) wide-area forecasting models.

[0734] 2. Dynamic Quantitative Reasoning

[0735] During the deployment phase of a neural network model, dynamic quantization can be used to convert model parameters from 32-bit floating-point numbers to 8-bit integers while maintaining high precision in the forward computation process. By quantizing linear and convolutional layers using `torch.quantization.quantize_dynamic`, the model size can be compressed to 1 / 4 of its original size, inference speed can be increased by 2-3 times, and accuracy loss can be kept within an acceptable range.

[0736] 3. Multi-resolution feature fusion

[0737] For forecasting tasks covering very large areas, a multi-resolution processing strategy can be adopted: first, the input is downsampled to extract global features at low resolution; then, it is fused with the original high-resolution features to balance computational efficiency and preservation of local details.

[0738] (vii) Online adaptive mechanism

[0739] To address the slow changes in the atmospheric environment and emission sources, this method supports an online adaptive update mechanism. Specifically, after model deployment, incremental learning can be achieved through a lightweight adaptation module:

[0740] Keep the backbone network parameters fixed and train only a lightweight adaptation layer (such as a low-rank adapter).

[0741] Maintain an experience replay buffer to store recent high-quality observation data;

[0742] Periodically fine-tune the adaptation layer using buffer data to enable the model to adapt to environmental changes while avoiding destructive forgetting.

[0743] This online adaptive mechanism enables the model to maintain prediction accuracy over long-term operation without the need for frequent retraining of the entire model.

[0744] The deterministic forecast mean field and covariance matrix generated by the method provided in this embodiment are used as wide-area background field data for probabilistic downscaling and uncertainty estimation of atmospheric environmental factor concentrations at local stations. For specific implementation, please refer to Embodiment Six.

[0745] Furthermore, the method provided in this embodiment can be widely applied to the following scenarios:

[0746] Real-time air quality forecasting and source analysis at the city level: providing environmental protection departments with high-precision, confidence-based forecasting products to support pollution source tracing and emergency decision-making;

[0747] Industrial park pollution diffusion risk assessment: rapid probability prediction of sudden leakage events and assessment of diffusion risk under different meteorological conditions;

[0748] Early warning of cross-border pollutant transport: predicting the trajectory and assessing the impact of long-distance transport events such as sandstorms and forest fire smoke;

[0749] Air quality simulation in carbon neutrality pathway planning: providing simulation tools for environmental impact assessment of carbon neutrality policies.

[0750] In terms of business model, this solution can provide regulatory tools to government environmental protection departments, route planning services to new energy vehicle companies, and personalized air quality early warning applications to the public, and has broad market application prospects.

[0751] The forecasting method provided in this embodiment differs from traditional deterministic forecasting methods in the following ways. Deterministic methods treat atmospheric forecasting as a function approximation problem, directly outputting a single numerical value in Euclidean space, without including information on the reliability of the prediction results. This method reconstructs the forecasting problem as conditional density estimation, outputting multivariate Gaussian distribution parameters defined in a probability measure space, so that the forecast results simultaneously include the most probable concentration value and its complete probabilistic information. This probabilistic modeling approach has the following technical features:

[0752] 1. Quantification of Uncertainty

[0753] By explicitly modeling the covariance matrix, the prediction uncertainty for each spatial location and each pollutant can be quantified, providing derived indicators such as confidence intervals and exceedance probabilities for risk decision-making. For example, the specific probability that PM2.5 concentration in a certain area will exceed the national secondary standard at a future time can be obtained.

[0754] 2. Spatial correlation perception

[0755] The off-diagonal elements of the covariance matrix can capture the correlation of prediction errors between different spatial locations, enabling the model to perceive regional forecast biases. When the forecast for a certain region is systematically overestimated, the covariance matrix can reflect this spatial correlation, providing prior information on uncertainties for subsequent local corrections and helping to improve the correction effect of downstream tasks.

[0756] 3. Adaptive weighting

[0757] In the negative log-likelihood loss function, the accuracy matrix itself acts as a weighting agent, giving the model more attention to the prediction accuracy in high-confidence regions while reducing the impact of high-noise regions. This adaptive weighting mechanism can optimize parameters overall and improve the model's robustness in scenarios with uneven signal-to-noise ratio distributions.

[0758] The deterministic forecast mean field and covariance matrix generated by this method can be used as wide-area background field data, providing a macroscopic physical framework and uncertain prior information for subsequent local site refined probabilistic prediction (see Example 6), thereby realizing a prediction chain from wide-area to local and from deterministic to probabilistic.

[0759] Example 6

[0760] This embodiment relates to a site-level atmospheric environmental factor concentration probability prediction method, which belongs to the probability output link of the local model in the dual-expert fusion prediction architecture of this invention. It is used to output the probability distribution parameters of any prediction location based on the wide-area background field and site observations, and to provide risk indicators for environmental decision-making.

[0761] When making decisions, environmental management departments are not only concerned with the specific numerical values ​​of pollution concentrations, but also with the risk level of the predicted results, such as PM2.5 concentrations exceeding the national secondary standard (75 μg / m³). 3 What is the probability of ( )? However, most existing technologies provide deterministic point predictions and cannot quantify the uncertainty of the prediction. A few methods involving uncertainty have obvious drawbacks: Monte Carlo Dropout (MC Dropout) requires multiple forward propagations during inference, has high computational costs, and mainly estimates cognitive uncertainty rather than the inherent randomness of the data; Bayesian neural networks have extremely high training and inference complexity, making them difficult to deploy in practice; quantile regression can provide a prediction interval, but it cannot output a complete probability distribution, and multiple quantile predictions may overlap (quantile crossover problem), violating the basic axioms of probability.

[0762] In spatiotemporal forecasting tasks such as air quality forecasting, providing a single forecast value is often insufficient to meet the needs of practical applications. Due to inherent model errors, uncertainties in input data (such as weather forecast errors), and the complexity of the atmospheric system itself, forecast results inevitably exhibit a certain range of fluctuation. City managers are more concerned with "how great the risk is" than "what the only possible outcome will be" when making decisions. They need to know the specific probability that the concentration of atmospheric environmental factors will exceed the warning threshold in order to make scientific decisions on the timing and intensity of emergency response measures.

[0763] Meanwhile, practical operational systems typically begin with a wide-area, coarse-resolution background field, such as a 10km × 10km numerical weather prediction or AI-generated forecast field. However, urban management and the public are concerned with air quality and its risks at a smaller scale (1km or even street level). This involves the classic problem of statistical downscaling: given a coarse-grid forecast, how to provide a fine-resolution, uncertain pollutant prediction for any location outside the grid center (referred to as a "sub-pixel" in this paper). This problem is made even more challenging by the sparseness of monitoring stations and the coexistence of observation and model errors.

[0764] Based on this, this embodiment provides a site-level probability prediction method for atmospheric environmental factor concentrations. The core idea of ​​this method is to directly output the probability distribution parameters of the atmospheric environmental factor concentrations at the target location through a parameterized network, rather than a single predicted value; to adopt a maximum likelihood estimation framework and use negative log-likelihood as the loss function for end-to-end training, enabling the model to simultaneously learn the predicted mean and standard deviation; and to design a multimodal heterogeneous information fusion strategy, which deeply fuses information such as coarse-resolution meteorological fields, observations from surrounding sites, and static geographic attributes with the target sub-pixel location as a condition, thereby achieving probability prediction with uncertainty for any spatial location.

[0765] like Figure 6 As shown, the method includes the following steps:

[0766] S600: Obtain the geographic coordinates of any predicted location within the target area.

[0767] The predicted location is any spatial point within the target area, including but not limited to:

[0768] Monitoring site locations: The actual locations where air quality monitoring equipment has been deployed. Historical observation data is available for these locations and can be used as supervisory signals for model training.

[0769] Any location without monitoring data: Any point within the target area where no monitoring equipment is deployed, such as residential areas, schools, hospitals, traffic intersections, and other sensitive locations of public concern;

[0770] Regular grid points: The center point of each grid after dividing the target area into regular grids (such as 1km×1km or 500m×500m), used to generate a continuous spatial distribution field;

[0771] User-defined locations: Any location dynamically specified according to application needs, such as the location of a pollution incident or a hotspot for complaints.

[0772] One of the core advantages of this solution is that it can output probability predictions with uncertainty for any of the above locations (including locations without historical observation data), breaking through the limitation of existing technologies that can only make predictions for existing monitoring stations.

[0773] The geographic coordinates are represented using a unified coordinate system to ensure accurate spatial positioning and alignment of cross-modal information. Specifically, at least one of the following forms may be used:

[0774] Latitude and longitude coordinates: Geographic coordinates expressed in longitude and latitude (such as the WGS84 coordinate system), which are suitable for large-scale cross-regional forecasting and are easy to align with global or regional scale data such as meteorological reanalysis fields and satellite remote sensing data;

[0775] Projected plane coordinates: Coordinates after latitude and longitude are projected onto a plane rectangular coordinate system (such as UTM projection, Gauss-Kriging projection), which are suitable for high-precision prediction of local areas and facilitate distance calculation and local area clipping.

[0776] Grid index coordinates: row and column indices in a predefined regular grid system, suitable for regular gridded prediction tasks, and easy to connect with gridded processing modules such as convolutional neural networks;

[0777] Relative position encoding: The offset relative to a reference point (such as the city center or monitoring station) within the target area, suitable for specific application scenarios that require relative position information.

[0778] This embodiment preferably uses latitude and longitude coordinates because they are the most directly aligned with the spatial data of meteorological reanalysis fields and wide-area background fields. In specific implementation, the latitude and longitude coordinates can be normalized (e.g., scaled to the [-1,1] interval) to facilitate neural network training.

[0779] Depending on the prediction task, the geographic coordinates can be obtained in the following ways:

[0780] Preset grid batch generation: For tasks that require generating continuous spatial distribution fields, the spatial range (such as minimum longitude, maximum longitude, minimum latitude, and maximum latitude) and spatial resolution (such as 0.01°×0.01° or 1km×1km) of the target area are preset, and the coordinate sequence of all grid points is automatically generated through grid traversal;

[0781] User-specified location input: Receive a list of coordinates of one or more predicted locations specified by the user through a user interface, API interface, or file upload, such as key schools and hospitals specified by the environmental protection department;

[0782] Monitoring site coordinate reading: The coordinates of existing sites are read directly from the monitoring site information database and used for prediction of the monitoring site itself or as a benchmark for model verification;

[0783] Event-driven coordinate acquisition: In response to pollution incident alarms or public complaints, dynamically acquire the coordinates of the incident location for real-time prediction.

[0784] Before inputting geographic coordinates into the subsequent parameterized network, necessary preprocessing and encoding of the coordinates can be performed to enhance the model's ability to represent spatial locations:

[0785] Coordinate normalization: Maps the original latitude and longitude coordinates to a uniform numerical range (such as [-1,1] or [0,1]) to avoid training instability due to an excessively large range of coordinate values;

[0786] Location embedding: By using learnable embedding layers or location encoding functions (such as sinusoidal location encoding, Gaussian kernel encoding), low-dimensional coordinate inputs are mapped to a high-dimensional feature space, enabling the model to capture spatial location differences more precisely.

[0787] Multi-scale positional encoding: Simultaneously employing positional encoding at multiple scales (such as absolute positional encoding and relative positional encoding) enables the model to perceive both global spatial structure and capture local spatial relationships;

[0788] Integration with static attributes: Geographic coordinates are combined with pre-acquired static geographic attributes (such as elevation and land use type) to form a richer location representation.

[0789] The obtained geographic coordinates will play a crucial role in subsequent steps:

[0790] In wide-area background field data fusion, it serves as a spatial index for extracting local background information from raster data, guiding coordinate conditional convolution or attention mechanism sampling.

[0791] In the dissemination of information about surrounding sites, it is used as a query condition to filter neighboring sites, and as a query vector in the cross-attention mechanism to dynamically aggregate historical information of surrounding sites.

[0792] In static geographic attribute association, it is used as a spatial key to look up the corresponding static attribute data;

[0793] In the parameterized network input, the pre-constructed spatiotemporal context vector is used together as input, enabling the network to output adaptive probability distribution parameters based on location differences.

[0794] Through the above design, this step provides a precise spatial positioning foundation for the entire method, ensuring that all subsequent information fusion and probability prediction can be customized based on the target location, thereby achieving refined uncertainty prediction for any location.

[0795] S610, acquire multimodal input information related to the predicted location, the multimodal input information including at least wide-area background field data, historical observation data of surrounding monitoring stations, and static geographic attribute data.

[0796] This step aims to collect and preprocess multi-source heterogeneous data related to the predicted target location, providing a comprehensive information foundation for subsequent probability distribution parameter prediction. This data comes from different modalities, spatial resolutions, and update frequencies, and needs to be integrated through a unified acquisition and preprocessing mechanism.

[0797] (a) Wide-area background field data

[0798] The wide-area background field data originates from the wide-area atmospheric environmental factor concentration probabilistic prediction method described in Embodiment 2 or Embodiment 5 of this invention, and its output is a wide-area coarse-resolution background field covering the entire target area. Specifically, it includes the following two forms:

[0799] Deterministic forecast field: A predicted field of average concentrations of atmospheric environmental factors for future time steps, such as the concentration distribution of pollutants like PM2.5, PM10, and O3. This data is stored in a regular grid format, with each grid point containing a deterministic concentration prediction value, characterizing the pollution transport background and meteorological forcing conditions at the regional scale.

[0800] The forecast uncertainty information in the form of a covariance matrix is ​​generated by the probabilistic forecasting method for wide-area atmospheric environmental factor concentrations described in Example 5, and characterizes the correlation of prediction errors of the wide-area model at different spatial locations. This covariance matrix can serve as prior information for local uncertainty estimation, constraining the generation of subsequent probability distribution parameters. When the uncertainty of the wide-area model is high in a certain region, the standard deviation of the probability distribution of the local model increases accordingly, reflecting a decrease in prediction confidence.

[0801] After obtaining the above data, spatial alignment and localization processing are required based on the geographic coordinates of the target's predicted location. The specific methods will be described in detail in subsequent steps.

[0802] (ii) Historical observation data from surrounding monitoring stations

[0803] Monitoring stations provide accurate local atmospheric environmental factor concentration observations, serving as a core data source for correcting and refining the wide-area background field. For each target prediction location, historical observation data from surrounding monitoring stations must be obtained, specifically including:

[0804] Site filtering: Based on the geographic coordinates of the predicted target location, monitoring stations within a preset distance range are filtered out. Filtering strategies can employ a fixed radius threshold (e.g., all stations within 10 kilometers), a dynamic radius (adaptively adjusted based on station density), or the K-nearest neighbor method (selecting the N1 nearest stations). This filtering process retains only stations with a strong spatial correlation to the target location, avoiding the inclusion of irrelevant stations that are too far away in the calculation, thus improving information utilization efficiency.

[0805] Historical time-series data acquisition: For each selected monitoring station, acquire its historical atmospheric environmental factor concentration series over a past period, including hourly concentration values ​​of various pollutants such as PM2.5, PM10, and O3. The historical time step length can be set according to the needs of the forecasting task, such as the past 24 hours, 48 ​​hours, or 72 hours.

[0806] Temporal feature encoding: Feature encoding is performed on the historical time-series data of each site to extract feature vectors that can characterize the historical pollution patterns of the site. The encoder can use temporal models such as long short-term memory networks, gated recurrent units, or one-dimensional convolutional neural networks to compress the variable-length historical sequence into a fixed-dimensional feature representation, which facilitates subsequent fusion with information from the target location.

[0807] To further enhance the representation capability of historical information of monitoring stations, the dynamic spatial context encoding method described in Embodiment 3 of this invention can be used to enhance the aforementioned temporal features. Specifically, the historical time-series data of each monitoring station is input into the dynamic spatial context encoding module. This module generates a spatial context encoding sequence containing the spatiotemporal interaction relationship between stations through parameterized transformation and attention mechanism, serving as an enhanced feature of the station's historical information. Compared to directly using the original time-series data or independently encoded features, this enhanced feature can more comprehensively reflect the pollution transmission relationship and dynamic change patterns between stations (e.g., the impact of pollutant transport changing with wind direction), providing a higher quality information foundation for the subsequent cross-attention aggregation step (see S622), thereby improving the accuracy of probability prediction at any location and the rationality of uncertainty estimation.

[0808] (III) Static geographic attribute data

[0809] Static geographic attribute data describes the inherent characteristics of the target's predicted location and its surrounding environment. These characteristics are relatively stable and do not change rapidly over time, but they have a significant impact on the spatial distribution of atmospheric environmental factor concentrations. These mainly include:

[0810] Basic geographic information: Topographic features such as elevation, slope, and aspect affect the conditions for pollutant diffusion;

[0811] Land use type: such as industrial land, commercial land, residential land, farmland, forest land, water body, etc., the emission intensity of different land use types varies significantly;

[0812] Road network density: Reflects the intensity of traffic emissions, and can be calculated by weighting the length of roads within a certain radius or the road class.

[0813] Population density: a measure of the intensity of human activity and is related to emissions from residential sources;

[0814] Emission source information: Distribution of emission intensity from industrial point sources and area sources;

[0815] Sensor attributes: For the location of the monitoring station, this also includes equipment parameters such as sensor type and measurement accuracy.

[0816] These static attribute data are typically stored in raster or vector form and can be spatially queried and extracted based on the geographic coordinates of the target's predicted location, forming a feature vector bound to the location.

[0817] S620, the multimodal input information and the geographic coordinates of the predicted location are input into the parameterized network, and the parameterized network outputs a parameter set of the probability distribution of the atmospheric environmental factor concentration at the predicted location in the future time step.

[0818] This step is the core of probabilistic downscaling and uncertainty estimation. To achieve probabilistic distribution prediction, this invention adds a parameterized network to the output of existing local site prediction models (especially the pluggable time-series prediction module described in Example 4). This parameterized network takes multimodal fusion information as input and directly outputs the probability distribution parameters of the atmospheric environmental factor concentrations at the target location, thereby extending deterministic prediction to probabilistic prediction.

[0819] This step specifically includes the following sub-steps:

[0820] S621, using the geographic coordinates of the predicted location as a condition, extract local background features from the wide-area background field data to generate a location-sensitive background context.

[0821] Wide-area background field data is a coarse-resolution background field covering the entire target area, containing deterministic forecast fields and forecast uncertainty information in the form of covariance matrices (see Example 2 or Example 5). Since the predicted target location is usually not exactly located at a grid point, and the impact of different spatial locations on the target varies, it is necessary to extract local information related to the target location from the wide-area background field data. Specifically, one or a combination of the following methods can be used:

[0822] Coordinate-aware sampling: This method directly extracts the values ​​of the grid points closest to the target location. It is computationally efficient and suitable for real-time prediction scenarios requiring rapid response. For higher accuracy, bilinear interpolation can be used to perform a weighted average of the surrounding four grid points.

[0823] Local region cropping and coordinate-conditional convolution: Centered on the target location, a local region of a preset size (e.g., a 5×5 or 7×7 grid) is cropped from the wide-area background field, preserving the neighborhood spatial structure. To enable the convolution operation to perceive the precise location of the target, the target coordinates (α, β) are embedded and encoded, and then concatenated with the cropped local feature map along the channel dimension through a broadcast mechanism. This concatenation is then fed into a series of 2D convolutional layers for processing. These coordinate-conditional convolutional layers can learn location-sensitive local meteorological feature representations, enabling the model to distinguish the differences in impact when the target is located at the center, boundary, or corner of the grid.

[0824] Attention mechanism sampling: At a finer level, the coded coordinates of the target location are used as the query, and the features of each grid point within the clipping region are used as the key and value. The weight of the target location to each point within the clipping region is calculated through an attention mechanism, and then weighted aggregation is performed. This approach allows the model to dynamically focus on the meteorological points that have the greatest impact on the target location within a local area. For example, when the target location is downwind of a strong emission source, the attention weights will automatically shift in that direction.

[0825] Feature Map Generation and Deformable Convolution: The entire wide-area background field is processed through a series of convolutional layers to generate multi-channel, high-level feature maps. Then, using methods such as deformable convolution, the offset of the feature vector corresponding to the target location is learned and accurately sampled from the feature maps. This method can achieve more refined spatial feature extraction, and the sampling position can be adaptively adjusted according to task requirements.

[0826] Using the above method, a local background context vector, i.e., a location-sensitive background context, is generated that can characterize the wide-area background features around the target location, denoted as f. background This vector not only contains background concentration information at the target location, but also incorporates the spatial structure features of the surrounding neighborhood, providing macroscopic physical constraints for subsequent probabilistic predictions.

[0827] S622, using a spatial attention mechanism or graph neural network, historical observation data from surrounding monitoring stations are propagated to the predicted location to generate a context vector reflecting the influence of the stations.

[0828] The goal of this step is to effectively transfer historical observation data from surrounding monitoring stations to the target prediction location, enabling subsequent probabilistic predictions to fully utilize the constraints of local real-world observations. Since the monitoring stations are discretely distributed, an intelligent information propagation mechanism is required, allowing the parameterized network to dynamically filter and aggregate information from surrounding stations based on the current context of the target location. Specifically, this includes the following sub-steps:

[0829] Nearest neighbor site selection: Based on the geographic coordinates of the target's predicted location (usually represented by longitude α and latitude β), monitoring stations within a preset distance range are selected. Selection strategies can employ a fixed radius threshold (e.g., all stations within 10 kilometers), a dynamic radius (adaptively adjusted based on station density), or the K-nearest neighbor method (selecting the N nearest stations). This selection process retains only stations with a strong spatial correlation to the target location, avoiding the inclusion of irrelevant stations that are too far away in the calculation, thus improving information utilization efficiency.

[0830] Temporal Feature Encoding: For each selected monitoring station, its historical air quality time series over a past period is obtained, and features are extracted using an independent temporal encoder. The temporal encoder can employ a long short-term memory network, a gated recurrent unit, or a one-dimensional convolutional neural network to compress the variable-length historical sequence into a fixed-dimensional feature vector. This feature vector characterizes the station's historical pollution patterns and temporal evolution, serving as the basis for subsequent information aggregation.

[0831] Cross-attention dynamic aggregation: In order to achieve intelligent filtering and fusion of information from surrounding sites, this step adopts the cross-attention mechanism described in Example 3 S340. The geographic coordinate code / local background context of the target predicted location is used as the query vector, and the historical feature vectors of surrounding sites are used as key and value vectors. Attention weights are calculated and weighted aggregation is performed. Unlike Example 3, the query vector in this step is bound to the spatial features of the predicted location, so as to realize the dynamic propagation of site information to any predicted location.

[0832] Specifically, we first construct three types of vectors:

[0833] Query vector: Generated based on the predicted target location, it can be the geographic coordinate encoding of that location or the local background context extracted from the wide-area background field in the previous step. The query vector represents the "information demand" or "current state" of the target location.

[0834] Key vector and value vector: For each surrounding site, the key vector is obtained by linearly transforming (or not transforming) its historical feature vector. The key vector is used to calculate the matching degree with the query vector, while the value vector carries the information of the sites to be aggregated.

[0835] Then, the weight of each site is determined through attention calculation: the similarity between the query vector and the key vector of each site is compared to obtain an initial attention score; then all scores are normalized so that the sum of the weights of all sites is 1. This weight reflects the importance of each site in the context of the current target location—sites with strong spatial correlation to the target location and whose historical patterns are consistent with the current context will be given higher weights, and vice versa.

[0836] Finally, the value vectors of each station are weighted and summed according to the weights mentioned above to obtain a context vector reflecting the station's influence. This vector integrates information from all surrounding stations, but highlights the contribution of the station most valuable to the current target location.

[0837] This dynamic aggregation mechanism can adaptively adjust site weights: for example, when the target location is downwind of a site, the weight of that site automatically increases; when a site's historical pollution pattern differs significantly from the current wide-area background (e.g., the site is affected by localized sudden emissions and the background field does not reflect this), its weight may decrease accordingly to avoid introducing inconsistent information.

[0838] As an alternative or supplement to the cross-attention mechanism, graph neural networks can also be used for information propagation. Specifically, a local graph is constructed with the target location and surrounding stations as nodes, and the edge weights between nodes are dynamically calculated based on geographical distance, wind direction, and other factors. Through the graph attention layer, each node can aggregate information from its neighboring nodes, and the edge weights are learned autonomously during training. After information propagation through multiple layers of graph neural networks, the target location node can draw on historical experience and spatial dependency information from its neighboring stations to generate a context vector that also reflects the influence of the stations.

[0839] Both methods can propagate site information to any location, and those skilled in the art can choose the appropriate technical path based on the actual application scenario. Through the above steps, the final generated site influence context vector is closely related to the target location, providing rich local observation constraint information for subsequent probability distribution parameter prediction.

[0840] S623, the location-sensitive background context, the context vector reflecting the site's influence, and the static geographic attribute features are fused to form a unified spatiotemporal context vector, which is used as the input to the parameterized network.

[0841] After extracting the above two modalities of information, it is necessary to integrate them with the static geographic attribute features of the target predicted location to form a unified spatiotemporal context vector:

[0842] Static geographic attribute encoding: The static geographic attribute data acquired by S610 is encoded. For numerical attributes (such as elevation and road network density), normalization is performed directly before use; for categorical attributes (such as land use type), they are mapped to low-dimensional dense vectors through an embedding layer; the features of each attribute are concatenated or fused to form a static feature vector f. staticfstatic .

[0843] Multi-source feature concatenation: Encoding the geographic coordinates p of the target location and the local background context f background Site Impact Context f stations Static features fstatic The features are then concatenated to form a high-dimensional feature vector.

[0844] Feature fusion network: Optionally, the concatenated vectors can be input into a lightweight fusion network for further processing. The fusion network can employ a multilayer perceptron structure, containing one or more fully connected layers and nonlinear activation functions, to achieve cross-modal information interaction and dimensionality compression.

[0845] The design of converged networks aims to enable information from different sources to calibrate and enhance each other, such as enabling site observation information to locally correct the wide-area background and enabling static geographic attributes to adapt the meteorological field to the terrain.

[0846] The spatiotemporal context vector constructed through the above steps comprehensively characterizes the complex environmental background of the target prediction location in future time steps, including multi-dimensional information such as wide-area transmission background, local observation constraints, and geographical environmental features, providing rich input features for the subsequent prediction of probability distribution parameters.

[0847] S624, the spatiotemporal context vector is input into the parameterized network, and after forward propagation, the probability distribution parameter set is output.

[0848] The parameterized network takes the spatiotemporal context vector constructed by S623 as input, performs forward propagation through a series of network layers, and finally outputs a parameter set of the probability distribution of atmospheric environmental factor concentrations at the target location in future time steps.

[0849] Network backbone: The parameterized network can adopt a multilayer perceptron structure, containing multiple hidden layers and non-linear activation functions (such as ReLU, GELU). The number and dimension of hidden layers can be adjusted according to the task complexity, with a typical configuration of 2-4 layers and 128-512 neurons per layer.

[0850] Output layer design: The number of neurons in the output layer depends on the number of free parameters of the chosen probability distribution. Taking the most commonly used normal distribution as an example, the output layer needs two neurons, corresponding to the mean μ and the standard deviation σ, respectively. For the log-normal distribution, the output layer also outputs two parameters, namely the mean μ in the logarithmic domain. log and standard deviation σ logIf multiple pollutants need to be predicted simultaneously, the output layer dimension is expanded to 2×C2×C, where C is the number of atmospheric environmental factors and C2 is the probability distribution feature dimension. In this embodiment, C2=2, corresponding to the mean and standard deviation of each pollutant's probability distribution parameters. If it is necessary to further characterize higher-order features such as skewness and kurtosis of the probability distribution, C2 can be expanded to 3 / 4. Those skilled in the art can adjust it according to actual needs. The final output layer dimension is C2×C×H×W (H / W is the spatial grid height / width of the prediction area), ensuring that each spatial grid point, each pollutant, and each probability distribution parameter has a corresponding output value.

[0851] Positive definiteness constraint of standard deviation: Since the standard deviation must be positive, the neurons in the output layer corresponding to this parameter need to be processed by activation functions that ensure a positive output. Common methods include:

[0852] Exponential function: σ = exp(z) σ The output is always positive and increases monotonically with the input.

[0853] Softplus function: σ = log(1 + exp(z)) σ The output is positive and close to 0 when the input is a large negative number, resulting in better numerical stability.

[0854] ELU+1:σ=ELU(z σ +1, the output range is greater than 1, and it is suitable for scenarios where there is a priori knowledge of the lower limit of standard deviation.

[0855] Fusion of Covariance Matrix Priors: If the wide-area background field data contains forecast uncertainty information in the form of a covariance matrix, this information can be used as a prior input into the parameterized network. Specifically, the covariance matrix elements (such as variance values) corresponding to the target location can be used as additional features, concatenated with the spatiotemporal context vector, and then input into the network, enabling the local uncertainty estimation to perceive the confidence level of the wide-area forecast. When the uncertainty of the wide-area model is high in a certain region, the parameterized network can correspondingly increase the standard deviation of the output to reflect the decrease in prediction confidence.

[0856] Through the above design, the parameterized network outputs a set of probability distribution parameters θ=(μ, σ) (taking the normal distribution as an example) for each target prediction location. These parameters fully describe the probability distribution characteristics of atmospheric environmental factor concentrations at the location in future time steps.

[0857] In this invention, the parameterized network is trained end-to-end by minimizing the negative log-likelihood loss function, which is calculated based on the probability density between the probability distribution parameters predicted by the parameterized network and the actual observed values.

[0858] The training objective of parameterized networks is to ensure that the probability distributions they predict can explain the true observed data to the greatest extent possible; this is the core idea of ​​maximum likelihood estimation. In deep learning, this objective is achieved by minimizing the negative log-likelihood loss function.

[0859] This invention employs negative log-likelihood as the training loss function, enabling the parameterized network to simultaneously learn the mean and standard deviation parameters of the predicted probability distribution during training. Taking the normal distribution as an example, the loss function includes a log-standard deviation term and a standardized squared residual term. The log-standard deviation term acts as a regularization term, preventing the model from overestimating uncertainty. The standardized squared residual term dynamically adjusts the optimization weights of the mean based on the currently estimated variance. When the variance is large, the penalty for the mean error by this term decreases, allowing the model to have greater prediction bias in regions of high uncertainty.

[0860] Based on the distribution characteristics of atmospheric environmental factor concentrations, a normal or log-normal distribution can be chosen as the probability distribution form for the predicted output. If the atmospheric environmental factor concentrations exhibit a right-skewed distribution, a log-normal distribution can be used for modeling. In this case, the parameterized network predicts the mean and standard deviation of the log-transformed variables, and the loss function form is the same as that of the normal distribution; only the true input values ​​need to be replaced with the log-transformed values.

[0861] To further improve the spatial consistency of the probabilistic prediction field, a consistency penalty term for the predicted probability distributions of adjacent locations can be introduced into the loss function. This constrains the mean or standard deviation of adjacent locations to change smoothly in space, enabling the parameterized network to output a smooth uncertainty estimate in space. For example, for two spatially adjacent prediction locations j1 and j2, the mean or standard deviation of their prediction distributions can be constrained to change smoothly in space:

[0862] L smooth =∑ （j1，j2）∈G ||μ j1 -μ j2 || 2 +λ0×||σ j1 −σ j2 || 2 .

[0863] Where G is the set of all adjacent pairs, and λ0 is the hyperparameter that balances the weights of the standard deviation and mean terms.

[0864] The total loss function becomes L=L NLL +β0×L smooth .

[0865] Wherein, β0 is the overall balance coefficient of the spatial smoothness constraint, which is used to adjust the proportion of this constraint term in the total loss.

[0866] Optimization process: End-to-end training is performed using gradient descent optimizers (such as Adam and SGD). The gradient of the loss function with respect to all learnable parameters in the parameterized network is calculated using the backpropagation algorithm, and the parameters are adjusted synchronously according to the optimizer's update rules. The training process is iterative until the model converges or reaches a preset stopping condition (such as the validation set loss no longer decreasing).

[0867] Through the end-to-end training described above, the parameterized network not only learns to predict the expected values ​​of atmospheric environmental factor concentrations, but also learns the uncertainty of the predictions, enabling the model to output a complete probability distribution and providing comprehensive support for subsequent risk decisions.

[0868] S630, Based on the probability distribution parameter set, construct the complete probability density function of the atmospheric environmental factor concentration.

[0869] This step is the output stage of probabilistic prediction. Its goal is to transform the probability distribution parameters output by the parameterized network into an interpretable and applicable complete probability density function, providing a mathematical basis for subsequent risk decision-making.

[0870] (I) Construction of the probability density function

[0871] Based on the probability distribution parameter set output by the parameterized network and combined with the selected probability distribution type, a complete probability density function of atmospheric environmental factor concentrations at the target prediction location in future time steps is constructed.

[0872] If a normal distribution is used, the parameterized network outputs the mean and standard deviation parameters. Based on these two parameters, a probability density function of the normal distribution can be constructed. This function describes the relative probability of atmospheric environmental factor concentrations taking various possible values. Its shape is determined by the mean, which determines the center of the distribution, and by the standard deviation, which determines the degree of dispersion.

[0873] If a log-normal distribution is used, the parameterized network outputs the mean and standard deviation parameters of the variables after logarithmic transformation. In this case, the probability density function of the normal distribution in the logarithmic domain is first constructed, and then the probability density function of the atmospheric environmental factor concentration at the original scale is obtained through variable transformation. The log-normal distribution is suitable for scenarios where the atmospheric environmental factor concentration exhibits a right-skewed distribution, ensuring that the predicted values ​​are always positive and better fitting the distribution characteristics of the high-concentration tail.

[0874] (II) Mathematical expression of probability density function

[0875] The constructed probability density function fully characterizes the uncertainty of atmospheric environmental factor concentrations at the target location at future times. This function possesses the following properties:

[0876] The integral is 1 over the entire range of values, which conforms to the basic axioms of probability distribution;

[0877] It can calculate the cumulative probability within any concentration range, such as the probability that the concentration of an atmospheric environmental factor is below a certain threshold.

[0878] It can extract numerical features of a distribution, such as expected value, variance, and quantiles.

[0879] (III) The Application Basis of Probability Density Function

[0880] Based on the constructed probability density function, various decision support information can be further derived:

[0881] Confidence interval: Based on the probability density function, the prediction interval at a given confidence level is calculated. For example, a 90% confidence interval means that the concentration of atmospheric environmental factors has a 90% probability of falling within this interval.

[0882] Probability of exceeding the standard: For a given environmental quality standard threshold (such as a daily average PM2.5 concentration of 75 μg / m³), calculate the cumulative probability of atmospheric environmental factor concentrations exceeding the threshold, providing a quantitative risk indicator for early warning and emergency decision-making;

[0883] Quantile prediction: Extracting specific quantiles (such as median, quintile, etc.) from the probability density function to provide more robust prediction reference values.

[0884] (iv) Output format

[0885] For a single predicted location, the constructed probability density function can be directly used for subsequent risk analysis. For the entire predicted area, the probability density function parameters of all predicted locations can be organized to form a probability forecast field, supporting spatial risk zoning and visualization.

[0886] Through this step, the abstract parameters output by the parameterized network are transformed into density functions with clear probabilistic meaning, so that the prediction results not only include "what is most likely to happen", but also "how likely it is to happen", providing more comprehensive decision support information for environmental management departments.

[0887] This embodiment, together with the aforementioned Embodiments 3 and 4, constitutes a complete site spatiotemporal prediction system. The overall workflow and training mechanism of the system are as follows:

[0888] I. Overall System Workflow

[0889] (a) Data Input Layer

[0890] The system receives three types of input data:

[0891] Historical observation data: historical atmospheric environmental factor concentration sequences (such as PM2.5, PM10, O3, etc.) and local meteorological observation data (temperature, humidity, wind speed, wind direction) of each monitoring station.

[0892] Wide-area background field data: The grid values ​​corresponding to each location and their forecast uncertainty information in the form of covariance matrix extracted from the forecast field output by the wide-area model;

[0893] Auxiliary information includes station identification information, static geographic attribute data (such as elevation, land use type), time characteristics (hour, day of the week, holiday identifiers), and future weather forecast data.

[0894] (II) Spatial Coding Layer (Example 3)

[0895] The input data is first fed into the dynamic spatial context encoding module. This module independently processes the multidimensional input features of all monitoring stations at each discrete time step, generating context codes that characterize the dynamic spatial interactions between stations through parameterized transformations and attention mechanisms. This process is repeated for multiple historical time steps, and the codes generated at each time step are combined chronologically to form a spatially encoded feature sequence. This sequence fully preserves the dynamic spatial interaction information of each station at historical time steps, providing high-quality feature input for subsequent time series prediction.

[0896] (III) Time Series Prediction Layer (Example 4)

[0897] The spatially encoded feature sequences are input to a pluggable time-series prediction module. This module features standardized input / output interfaces and supports a hybrid expert architecture. In hybrid expert mode, multiple expert prediction subunits process the input sequences in parallel. A gating network dynamically generates expert weights based on site identification information and input sequence features. The outputs of each expert are weighted and fused to obtain the predicted concentrations of atmospheric environmental factors for each monitoring station at future time steps. This module can perform specialized modeling for different site types or different pollution processes, achieving adaptive prediction.

[0898] (iv) Uncertainty estimation layer (in this embodiment)

[0899] The parameterized network described in this embodiment is connected to the output of the time-series prediction module. For each prediction location (including monitoring stations and arbitrary sub-pixel locations), the network takes the following information as input:

[0900] Spatiotemporal context vectors obtained from the spatial coding layer;

[0901] The geographic coordinates of the target location;

[0902] Local background features and covariance matrix information extracted from a wide-area background field;

[0903] Static geographic attribute characteristics.

[0904] The parameterized network outputs a set of parameters (such as the mean and standard deviation of a normal distribution) for the probability distribution at each target location and each prediction time step through forward propagation. These parameters fully describe the probability distribution characteristics of future atmospheric environmental factor concentrations at that location.

[0905] (v) Output layer

[0906] The system ultimately outputs the complete probability distribution parameters for each predicted location, and can further construct probability density functions based on these parameters to calculate derivative risk indicators such as confidence intervals and out-of-range probabilities, forming a complete output chain from deterministic prediction to probabilistic quantification.

[0907] II. End-to-end training mechanism

[0908] The entire system is optimized using an end-to-end joint training approach. The training process is as follows:

[0909] (I) Loss Function Design

[0910] The total loss function consists of the following parts:

[0911] Negative log-likelihood loss: Calculated based on the probability distribution parameters output by the parameterized network and the actual observations at the monitoring stations, enabling the model to learn to predict both the mean and standard deviation simultaneously;

[0912] Spatial smoothness constraint: Introducing a consistency penalty term for the probability distribution between adjacent prediction locations to ensure that the probability prediction field changes smoothly in space and improves physical rationality.

[0913] (II) Forward Propagation and Loss Calculation

[0914] For each training batch, the system sequentially performs forward propagation of spatial encoding, temporal prediction, and uncertainty estimation to obtain the probability distribution parameters of each monitoring station; calculates the negative log-likelihood loss based on the real observations, and calculates the spatial smoothness constraint term as needed; and obtains the total loss by weighted summation.

[0915] (III) Backpropagation and parameter update

[0916] The gradient of the total loss with respect to all learnable parameters in the system is calculated using the backpropagation algorithm, including:

[0917] The parameters of the first, second, and third transforms in the spatial coding module;

[0918] The parameters of each expert prediction subunit and the parameters of the gating network in the time series prediction module;

[0919] Parameters of a parameterized network.

[0920] A gradient descent-type optimizer is used to update all parameters synchronously, iterating until the model converges.

[0921] Through this end-to-end joint training, the spatial coding module, the temporal prediction module, and the uncertainty estimation module can be optimized in synergy: the spatial coding module learns to generate spatial features that are beneficial to temporal prediction and uncertainty estimation, the temporal prediction module improves prediction accuracy with the support of spatial features, and the uncertainty estimation module learns the expected value and volatility of the prediction from the data at the same time, ultimately forming a complete prediction system from macro background to micro site, and from deterministic prediction to probabilistic quantification.

[0922] The uncertainty estimation scheme introduced in this embodiment brings the following technical features to the site spatiotemporal prediction system:

[0923] (1) Provide probability distribution output to enrich decision-making information

[0924] The forecast output is expanded from a single numerical value to a complete probability distribution. Decision-makers can obtain the predicted values ​​and their corresponding probability information, such as the specific probability that the concentration of atmospheric environmental factors will exceed the warning threshold, providing more complete data support for decision-making.

[0925] (2) Support risk management and early warning functions

[0926] It can generate spatial distributions of prediction intervals and probabilities of exceeding standards. In scenarios such as air pollution emergency response, resource allocation, and public health early warning, it can be used to identify risk areas with a high probability of heavy pollution in the future, assisting in early warning and resource allocation.

[0927] (3) Enhance model interpretability

[0928] Analyzing probability distribution parameters (especially standard deviation) can reveal the predictive confidence level of the model under different spatiotemporal scenarios. A higher standard deviation may reflect insufficient model understanding of the scenario, high uncertainty in the input data, or significant environmental fluctuations, providing clues for model improvement or human intervention.

[0929] (4) Improve the robustness and generalization ability of the model

[0930] Probabilistic modeling forces the model to learn the inherent variability of the data, enabling it to form relatively robust feature representations when faced with highly uncertain inputs in the real world. Simultaneously, this approach helps the model generalize better to areas with sparse or unobserved monitoring stations.

[0931] The multimodal fusion wide-area atmospheric environmental factor prediction model (hereinafter referred to as the wide-area model) described in Embodiment 2, together with the local site spatiotemporal prediction model (hereinafter referred to as the local model) described in Embodiments 3, 4 and 6, constitutes an atmospheric environmental factor prediction system, which is used to realize multi-level prediction of the target area from the wide-area scale to the local site.

[0932] Specifically, the wide-area model takes field-space data or point-field fusion features after fusion processing as input, and generates a wide-area coarse-resolution background field covering the entire target area through a three-dimensional spatiotemporal sensing encoder and decoder. It can also optionally output the covariance matrix to quantify forecast uncertainty. This background forecast field reflects the transport, diffusion, and evolution of pollutants at the regional scale, providing a macroscopic physical framework and meteorological forcing background for local site forecasts.

[0933] The local model takes historical time-series data from monitoring stations as its core, integrates background field information and its covariance matrix output from the wide-area model, captures time-varying spatial correlations between stations through a dynamic spatial context encoding module, and then uses a hybrid expert time-series prediction module to perform specialized modeling of pollution patterns at different stations. Finally, it outputs the concentration prediction values ​​for future times at each station, and can be extended to probabilistic forecasts for any spatial location.

[0934] Through the aforementioned collaborative design, the wide-area model and the local model form a fusion prediction architecture of "macro-scale background transmission + micro-scale fine correction": the background field output by the wide-area model provides regional-scale evolution trend constraints for the local model, while the local model uses site observations to perform local corrections on the background field and feeds the correction information back to the utilization process of the wide-area background through an attention mechanism. Simultaneously, the uncertainty quantification results of the wide-area model provide prior information for local probability prediction, realizing the transmission and fusion of uncertainty in cross-scale prediction.

[0935] This prediction system combines the ability of a wide-area model to characterize regional transport processes with the ability of a local model to capture local features. It differs from single-model or traditional numerical simulation methods in terms of computational efficiency, prediction accuracy, and richness of decision support information. It can be applied to scenarios such as urban air quality forecasting, pollution process early warning, and regional transport contribution analysis.

[0936] This invention also provides an electronic device, including: at least one processor; and a memory communicatively connected to the at least one processor; wherein the memory stores instructions executable by the at least one processor, the instructions being configured to perform the method described in this invention.

[0937] This invention also provides a computer-readable storage medium storing computer-executable instructions for performing the methods described in this invention.

[0938] It should be understood that the various forms of processes shown above can be used to reorder, add, or delete steps. For example, the steps described in this invention can be executed in parallel, sequentially, or in different orders, as long as the desired result of the technical solution disclosed in this invention can be achieved, and this is not limited herein.

[0939] The specific embodiments described above do not constitute a limitation on the scope of protection of this invention. Those skilled in the art should understand that various modifications, combinations, sub-combinations, and substitutions can be made according to design requirements and other factors. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of this invention should be included within the scope of protection of this invention.

Claims

1. A method for probabilistic prediction of wide-area atmospheric pollutant concentration, characterized by, The method is implemented through a pre-built neural network model, which includes a three-dimensional spatiotemporal attention encoder, a mean prediction head, and a covariance prediction head. The method includes the following steps: Obtain the input spatiotemporal cube of historical time steps, wherein the input spatiotemporal cube contains multiple time steps, multiple spatial locations and multiple atmospheric parameters; The input spatiotemporal cube is fed into an encoder based on a three-dimensional spatiotemporal attention mechanism to extract spatiotemporal feature maps; Based on the spatiotemporal feature map, a deterministic forecast mean field for future times is generated through a mean prediction head; Based on the spatiotemporal feature map, a covariance matrix for forecasting uncertainty is generated through the covariance prediction head; The deterministic forecast mean field and the covariance matrix together define the multivariate Gaussian distribution of the atmospheric state at future times.

2. The method of claim 1, wherein, The encoder based on the three-dimensional spatiotemporal attention mechanism adopts a three-dimensional window attention mechanism, which expands the two-dimensional spatial window into a three-dimensional window that includes the time dimension. Self-attention calculation is performed within the non-overlapping three-dimensional window, while capturing the dynamic correlation of the atmosphere in the horizontal space, vertical layering and time dimension.

3. The method of claim 2, wherein, The dimensions of the three-dimensional window include the size of the time window, the size of the spatial height window, and the size of the spatial width window. By sliding the three-dimensional window on feature maps at different levels, hierarchical extraction of multi-scale spatiotemporal features can be achieved.

4. The method of claim 1, wherein, The covariance matrix is ​​modeled using a constraint method that combines block diagonalization and low-rank global correlation: The diagonal section is used to model the local correlations between different atmospheric parameters within each spatial location; The low-rank global correlation component is used to model long-range correlations across spatial locations. It compresses the high-dimensional covariance matrix into a product of low-rank factors through low-rank decomposition. An additional diagonal perturbation term is used to ensure the positive definiteness of the covariance matrix.

5. The method according to claim 4, characterized in that, The low-rank global correlation component is achieved by multiplying a learnable low-rank matrix by its transpose, where the rank of the low-rank matrix is ​​much smaller than the total number of spatial locations.

6. The method according to claim 1, characterized in that, The method further includes a physical constraint enhancement optimization step: during the parameter optimization process of the neural network model, a mass conservation constraint term is introduced, and the prediction results that do not conform to physical laws are penalized by calculating the divergence of pollutant flux.

7. The method according to claim 6, characterized in that, The mass conservation constraint term is calculated based on the predicted mean field and the input wind field data, and is used to guide the prediction results to follow the physical laws of atmospheric pollutant transport during the training process.

8. The method according to claim 1, characterized in that, The method employs a progressive multi-task training framework to optimize the neural network model in stages: the first stage trains only the mean prediction head to optimize the deterministic prediction accuracy; the second stage freezes the mean prediction head and trains the covariance prediction head to learn uncertainty modeling; the third stage jointly fine-tunes all modules and adds auxiliary prediction tasks when necessary.

9. An electronic device, characterized in that, Including processor and memory; The processor executes the steps of the method as described in any one of claims 1 to 8 by invoking programs or instructions stored in the memory.

10. A computer-readable storage medium, characterized in that, The computer-readable storage medium is used to store a program or instructions that cause a computer to perform the steps of the method as described in any one of claims 1 to 8.

Images