An urban carbon flux inversion optimization method based on joint scheduling and matrix dimension reduction

Abstract

The present application relates to the field of atmospheric environment monitoring, and specifically discloses a kind of urban carbon flux inversion optimization method based on joint scheduling and matrix dimensionality reduction, comprising: ensemble sample generation scheduling step, through serial-parallel hybrid scheduling strategy, continuously generate multiple samples in single model process, reduce model call redundancy overhead;Covariance matrix dimensionality reduction calculation step, use Kronecker product to decompose prior error covariance matrix into space-time component, combined with multi-thread parallel architecture to build efficient computing link;And carbon flux inversion solving step, based on Bayesian set integration framework, the foregoing optimization result is fused to complete inversion.The present application significantly reduces the consumption of computing resources and memory occupation of urban high-resolution carbon flux inversion through the collaborative optimization of calculation scheduling mechanism and storage architecture, can realize the efficient and stable inversion of urban scale kilometer level, hour level carbon flux under general regular computing environment, effectively improves the real-time performance and engineering landing feasibility of carbon monitoring business.

Description

Technical Field

[0001] This invention belongs to the field of atmospheric environment monitoring, specifically a method for urban carbon flux inversion and optimization based on joint scheduling and matrix dimensionality reduction. Background Technology

[0002] Urban greenhouse gas carbon flux inversion is a key technical means to assess urban carbon emissions and support the achievement of carbon neutrality goals. This method couples atmospheric transport models with observational data to invert the spatiotemporal distribution of surface carbon flux based on a Bayesian ensemble assimilation framework. Among these methods, Monte Carlo assimilation methods such as ensemble Kalman filtering have become one of the mainstream methods for carbon flux inversion at regional and urban scales due to their ability to effectively handle nonlinear and non-Gaussian problems.

[0003] In this type of method, a certain size of ensemble sample needs to be generated through an atmospheric transport model to construct the observation operator and the error covariance matrix. As urban observation networks become increasingly dense and the spatiotemporal resolution of models continues to improve, the dimensionality of the state variables in the inversion system increases dramatically, and the size of the covariance matrix grows quadratically, causing storage requirements to rapidly exceed the memory capacity of conventional servers. At the same time, high-resolution atmospheric transport models themselves are computationally expensive, and traditional parallel strategies require independent model calculations for each ensemble sample, resulting in redundant runs of basic modules such as meteorological fields, significantly wasting computational resources and severely restricting the real-time performance and scalability of urban carbon flux inversion systems in operational scenarios.

[0004] Therefore, how to significantly reduce the computational redundancy of ensemble sample generation and the storage overhead of the covariance matrix while maintaining inversion accuracy has become a key technical problem that urgently needs to be solved to achieve real-time and operational inversion of high-resolution carbon flux in cities. Summary of the Invention

[0005] To address the aforementioned technical problems, this invention provides an optimization method for urban carbon flux inversion based on joint scheduling and matrix dimensionality reduction. This method solves the problems of high computational resource consumption, high storage requirements, and difficulty in achieving real-time operational operation caused by redundancy in ensemble sample generation and covariance matrix operation in existing high-resolution inversion systems.

[0006] This invention provides an optimization method for urban carbon flux inversion based on joint scheduling and matrix dimensionality reduction, comprising the following steps: Step S1, Collection Sample Generation Scheduling: Determine the total number of samples in the set. ,Will Configured for serial number generation With the number of parallel calls The product of, i.e. ; In a single atmospheric transport model invocation process, the code of the atmospheric transport model is modified to embed set dimension parameters, and a serial loop mechanism is used to execute only the target module in the model. In each iteration, an independent perturbation input parameter is generated. Concentration simulation samples are generated in a single model process call. One sample; Simultaneously, the job is submitted and executed concurrently through the supercomputing job scheduling system. Each of the aforementioned model invocation processes is independent; Step S2, Covariance matrix dimensionality reduction calculation: Based on the assumption of spatiotemporal separability, the prior error covariance matrix of the state variables is... Decomposed into time covariance matrix With spatial covariance matrix The Kronecker product, i.e. ,in, This represents the Kronecker product operator. The dimension is , The dimension is , For time steps, The number of spatial grid cells; The observation operator matrix Divide into blocks according to the time dimension, represented as follows: Each sub-block The dimension is , The total number of observations; In a shared memory architecture, a multi-threaded parallel computing approach is used, which involves traversing the matrix. For each column, calculate the corresponding matrix product. The column blocks are used directly in subsequent inversion solutions, thus avoiding the need to explicitly store the complete... A matrix, where each thread is responsible for calculation. Zhongyu A corresponding time column block; Step S3, Carbon flux inversion solution: Based on the Bayesian ensemble inversion framework, a cost function containing observation and background terms is constructed, and the observation operator matrix is ​​calculated using the ensemble samples generated in step S1. The linearized expression, and the result obtained using step S2. The matrix is ​​used to solve for the posterior values ​​of the state variables that minimize the cost function, thus completing the high-resolution inversion of urban carbon flux.

[0007] Preferably, the serially generated number The value range is greater than 1 and less than or equal to 50.

[0008] Preferably, the serially generated number The value of is determined through preliminary experiments and is used to balance the running time and computational resource consumption of a single model process. Its optimized value is 15 to 30.

[0009] Preferably, the number of parallel calls The value is dynamically determined based on the total number of available computing nodes and the computing resources of each node, ensuring that the total number of samples in the set is within a certain range. It meets the inversion accuracy requirements.

[0010] Preferably, the atmospheric transport model is the WRF-Chem model or a derivative thereof, and the modification of the model code includes adding controls. Input parameters and internal loop structure for generating concentration field set samples.

[0011] Preferably, the supercomputing job scheduling system is one of Slurm, PBS, or LSF, and the concurrent execution is achieved by combining the batch job script of the selected scheduling system with the MPI parallel environment.

[0012] Preferably, the multi-threaded parallel computing is implemented using OpenMP, POSIX Threads, or a multi-threading library built into the programming language.

[0013] Preferably, the state variable is urban surface carbon flux, and its dimension is... The observation data was acquired from ground-based and / or satellite sources. Concentration data, of which ground-based observation data are station-based discrete observations with a temporal resolution better than or equal to 1 hour and a spatial resolution better than or equal to kilometers.

[0014] Preferably, in step S2, the matrix operation task is assigned to the GPU accelerator for execution by extending the computing interface, so as to further improve computing efficiency.

[0015] Compared with the prior art, the present invention has the following beneficial effects: Through the collaborative innovation of two core technologies—serial-parallel joint scheduling and matrix Kronecker product dimensionality reduction parallel computing—the computational paradigm of traditional carbon assimilation inversion has been reconstructed at the system level. At the computational level, this invention breaks through the rigid parallel mode of "one model call corresponds to one sample" in traditional ensemble sample generation and innovatively designs a serial-parallel hybrid scheduling architecture. By embedding ensemble dimension control and cyclic generation mechanisms into the atmospheric transport model, it realizes the continuous generation of multiple perturbation samples in a single model process, fundamentally eliminating the redundant calculation of basic modules such as meteorological fields. Combined with the intelligent job distribution of the external supercomputing system, a hybrid execution flow of "internal serial and external parallel" is formed, which significantly reduces the number of model calls and the occupation of computing nodes while ensuring the statistical characteristics of the ensemble, and significantly alleviates the pressure of high-performance computing resources.

[0016] At the storage and computation acceleration level, this invention addresses the bottleneck of covariance matrix operations caused by the explosion in scale during Bayesian set homogenization by introducing a dimensionality reduction strategy based on the Kronecker product structure, which is based on the assumption of spatiotemporal separability. By decomposing the high-dimensional covariance matrix into a tensor product of time and spatial components, the explicit storage and manipulation of ultra-large-scale dense matrices are avoided. Furthermore, by combining shared memory multi-threaded parallel technology, the dimensionality-reduced matrix operation task is efficiently mapped to a multi-core CPU architecture, achieving a dual reduction in computational complexity and memory usage. Attached Figure Description

[0017] Figure 1 This is a diagram illustrating the overall technical aspects of the computation-storage optimization for the urban greenhouse gas carbon assimilation inversion model of this invention. Figure 2 This is a schematic diagram of serial-parallel scheduling using WRF-Chem as the atmospheric transport model in an embodiment of the present invention; Figure 3 This is a diagram illustrating the serial-parallel scheduling efficiency using WRF-Chem as an example in this embodiment of the invention. Figure 4 This is a time-series flowchart for the multi-period rolling optimization of urban carbon assimilation in this invention; Figure 5 This is a comparative analysis diagram of computation and storage in a certain city's assimilation and inversion system, using Kronecker product dimensionality reduction in this embodiment of the invention. Figure 6 This is an information table of nine high-precision monitoring stations in 2024, as described in this embodiment of the invention. Detailed Implementation

[0018] 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.

[0019] Figure 1This illustrates the overall technical framework for the computation-storage optimization of the urban greenhouse gas carbon assimilation inversion model of this invention. For example... Figure 1 As shown, the core of this invention lies in the collaborative work of a joint scheduling module for serial and parallel computing and a module for matrix Kronecker product dimensionality reduction and parallel computing. The serial and parallel scheduling module is responsible for efficiently generating ensemble simulation samples of the atmospheric transport model, and its output serves as the input to the observation operator. The matrix dimensionality reduction and parallel computing module optimizes the large-scale covariance matrix operations in the Bayesian ensemble inversion framework. The outputs of both modules are jointly input to the carbon flux inversion solution module, ultimately obtaining high spatiotemporal resolution urban carbon flux inversion results. This framework comprehensively covers the entire process from data input and computational acceleration to result output, systematically solving the computational and storage bottlenecks in high-resolution inversion.

[0020] like Figure 4 As shown, the Figure 1 The process is shown for three assimilation cycles (cycle 1, cycle 2, and cycle 3). Each assimilation window (containing a 1-day assimilation window and a 4-day lag window) has a 5-day state variable, typically composed of the lag window and the assimilation window. The green-framed state variable represents the prior flux, the orange-framed represents the observed data, and the blue-framed represents the posterior flux after assimilation from the observed data. The flux symbol xi represents the flux set at time i, and the number in parentheses indicates how many times the flux has been optimized. The blue dots represent the corresponding time intervals. Concentration. The assimilation system operates as follows: 1) Using WRF-Chem, combined with the initial field The prior flux within period 1, running S sets Transmission, until Deadline; 2) Introduce the observed values ​​y(t+5) within period 1 into the assimilation system to assimilate the prior; 3) Using WRF-Chem, the posterior at time t within period 1 is used as the state variable to obtain the posterior at the next time (t+1). Initial concentration field; 4) The posterior flux from (t+1) to (t+5) in period 1 is used as the prior flux in period 2. The prior flux at (t+6) in period 2 is generally obtained by the dynamic model or provided by external flux product data. 5) The set of prior fluxes at time (t+6) in period 2 is generated by the given background error covariance, while the set of prior fluxes from time (t+1) to (t+5) comes directly from period 1. 6) Introduce the observation y(t+6) at time t+6, and assimilate the prior flux of period 2 to obtain the posterior flux. The above process is repeated for period 3 until the end of the study period.

[0021] Example: Carbon flux inversion in a certain city in 2024 This embodiment selects a city as the study area, within which nine high-precision [databases / locations] are set up. Concentration monitoring stations; site information can be found at [link to site information]. Figure 6 The observation data has a time resolution of 1 hour; this embodiment uses one assimilation period as an example; according to Figure 4 By executing the algorithm on a rolling basis over multiple periods, the daily average urban carbon flux of a city in 2024 with a spatial resolution at the kilometer level can be obtained. In this embodiment, the number of spatial grids is 111×81.

[0022] Step S1: Set Sample Generation Scheduling The goal of this step is to efficiently generate the total number of samples for the inversion ensemble. Traditional methods require calling Sub-atmospheric transport models lead to a large amount of redundant meteorological field calculations; this invention adopts a joint serial-parallel scheduling strategy.

[0023] First, determine the total number of samples in the set. ;Will Split into serially generated numbers and number of parallel calls ,satisfy .

[0024] Serial generation implementation: Within a single atmospheric transport model call process, the model code is modified to embed set dimension parameters, and the model is run continuously using a loop mechanism. Each iteration perturbs the input parameters, thereby generating [the necessary parameters] during a single model process startup. An independent Concentration simulation samples; this process eliminates the overhead of redundant meteorological field calculations.

[0025] Parameter optimization and determination: The value of the serial number s1 was determined through preliminary experiments to balance time efficiency and resource utilization; such as Figure 3 As shown, within a research area of ​​a certain city, 28 computing cores were used to call the WRFChem model in a single run, and tests were conducted to generate different quantities of data. The total time consumed by the concentration field set and the average time consumed by a single set were compared. Experimental results show that the average generation time of a single set decreases power-lawfully with the increase of the number of sets. This invention quantitatively determines the gradual flattening of the decay trend by the relative decrease in the average time consumed by a single set corresponding to an adjacent number of sets. The formula for calculating the relative decrease is as follows:

[0026] In the formula, For the first Average running time of a single set corresponding to a number of serial numbers; The average time consumed by a single set when the serial number is i+1 is given; when the relative decrease R of two consecutive adjacent sampling points is less than the preset threshold (15% for example), the decay trend is determined to be flat.

[0027] Based on the results of this preliminary experiment, when the number of ensembles is 15, 20, 25, and 30, the average time per ensemble is 14.27 min, 10.67 min, 9.47 min, and 8.20 min, respectively. Increasing the number of ensembles from 15 to 20 results in a relative decrease of approximately 25.2%; from 20 to 25, approximately 11.2%; and from 25 to 30, approximately 13.4%. The relative decreases in the two consecutive groups are all less than the 15% threshold, indicating a gradual decline in performance. Therefore, this embodiment selects s1=20 as the optimal value, which minimizes the number of model starts within an acceptable time cost. This preset threshold can be adaptively adjusted according to the number of computational cores and the simulation study area, and is not limited to 15%.

[0028] Parallel scheduling implementation: based on Calculate the number of parallel calls Three independent jobs are submitted concurrently through the supercomputing job scheduling system. Each job executes a modified WRF-Chem model process, and each process serially generates 20 samples. The supercomputing scheduling system is the Slurm system, which combines an MPI parallel environment to achieve concurrent job execution. The scheduling process is illustrated below. Figure 2 As shown.

[0029] Through the above scheduling, a task that traditionally required 60 independent parallel model calls is optimized to require only 3 parallel calls. For example... Figure 3 As shown, this strategy significantly reduces the number of computing nodes and model startup overhead while maintaining accuracy. Combined with subsequent matrix dimensionality reduction techniques, such as... Figure 5 As shown, the computational and storage efficiency of the entire inversion system has been significantly improved.

[0030] Step S2: Covariance Matrix Dimensionality Reduction and Parallel Computation This step focuses on the covariance matrix, a core component in Bayesian ensemble inversion. The assimilation window and related calculations are optimized. In this embodiment, the assimilation window is 1 day, but the assimilation period contains a 5-day lag window (the lag window refers to the flux influence time window set before the observation time to ensure that the CO2 concentration at the observation time fully responds to the carbon flux over a previous period. Since atmospheric transport takes a certain amount of time, flux lags from multiple time periods are generally needed to jointly determine the concentration at a specified time point). (Concentration), the state variable is urban surface carbon flux, its total:

[0031] in, For time steps, Let P be the number of spatial grids; then one assimilation window P( )and ( Storage requirements are approximately 15GB and 370MB respectively, which is:

[0032]

[0033] One element is stored in 8 bytes; byte, .

[0034] In this embodiment, the assimilation period is 5 days, and there are 9 monitoring stations in the study area. The time resolution of the observation data is 1 hour, that is, 24 observations per day. Therefore, the total number of observations is n = 5 × 9 × 24 = 1080. HP storage is as follows: Figure 5 As shown in (a) “Without using Kronecker product”, the 1-year P and HP storage far exceeds the capacity of a regular server (approximately ~5.4T and ~132G respectively).

[0035] Matrix Dimensionality Reduction: Based on the assumption of spatiotemporal separability, the prior error covariance matrix is ​​reduced. Decomposed into time covariance matrix With spatial covariance matrix The Kronecker product, i.e.:

[0036] in, It is a 5×5 matrix. It is a (111×81)×(111×81) matrix; observation operator matrix The dimension is ,in To represent the total number of observations, the assimilation period in this embodiment is 5 days, with 9 monitoring stations in the study area and a time resolution of 1 hour for the observation data. ;Will The time-dimension representation is as follows Each sub-block The dimension is .

[0037] parallel computing Matrix: In shared memory, a multi-threaded parallel computing approach is used to compute the time covariance matrix. Each column is assigned to a separate computation thread; each thread is calculated based on its assigned value. Column vectors, spatial covariance matrix and the observation operator matrix The corresponding blocks are calculated. The corresponding column block in the matrix is ​​used directly for subsequent inversion solutions, thus avoiding the need to explicitly store the complete matrix. matrix.

[0038] This multi-threaded parallel computing is implemented using OpenMP, with the number of threads matching the number of processor cores; through the aforementioned dimensionality reduction and parallelization strategies, The computational complexity of matrices has decreased from traditional Significantly reduced, and further accelerated using multi-core parallelism.

[0039] like Figure 5 (b) and Figure 5 As shown in (c), after adopting this method, the number of assimilation windows within each window... The computation time was reduced from 94 seconds using the traditional method to 0.01 seconds, and the efficiency of posterior covariance calculation was also improved by orders of magnitude. Since explicitly storing the complete H·P matrix is ​​avoided, only the dimension-reduced D, E, and a small number of column blocks need to be stored, significantly reducing HP-related storage requirements (e.g., ...). Figure 5 (a)), fully compatible with regular servers.

[0040] Step S3: Solving for carbon flux inversion Based on the set of samples generated in step S1, the observation operator is calculated. The linearized expression; combined with the efficient calculation obtained in step S2. Substitute the matrix into the Bayesian set inversion framework.

[0041] Within the Bayesian ensemble inversion framework, a cost function is constructed that includes both observation and background terms:

[0042] in, For observation Concentration vector The observation error covariance matrix, The background carbon flux vector.

[0043] By minimizing the above cost function, the posterior estimate of the state variable can be obtained. Since step S2 has efficiently processed the data containing... The core matrix operations are simplified, thus greatly reducing the burden of solving this minimization problem.

[0044] Ultimately, the daily average urban carbon flux of a certain city for the entire year of 2024 was obtained with a spatial resolution of kilometers.

[0045] Effects and scalability analysis: like Figure 5 As shown in (d), the technical solution of the present invention has excellent scalability in terms of computational complexity and storage requirements compared with traditional methods; the performance advantage of the present method is more significant when the problem size increases.

[0046] In addition, to further improve computational efficiency, this method also includes an extended computing interface, which is used to allocate the matrix operation task in step S2 to the GPU accelerator for execution. By utilizing the parallel computing capabilities of the GPU, the inversion time can be further shortened.

[0047] This invention, through the organic combination of serial-parallel joint scheduling and matrix Kronecker product dimensionality reduction parallel computation, can achieve efficient and high-precision inversion of city-scale kilometer-level and hourly carbon flux on ordinary computing servers, providing a feasible technical solution for real-time monitoring and scientific assessment of urban carbon emissions.

[0048] The above description is merely a preferred embodiment of the present invention and is not intended to limit the scope of protection of the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A method for optimizing urban carbon flux inversion based on joint scheduling and matrix dimensionality reduction, characterized in that, Includes the following steps: Step S1, Collection Sample Generation Scheduling: Determine the total number of samples in the set. ,Will Configured for serial number generation With the number of parallel calls The product of, i.e. ; In a single atmospheric transport model invocation process, the code of the atmospheric transport model is modified to embed set dimension parameters, and a serial loop mechanism is used to execute only the target module in the model. In each iteration, an independent perturbation input parameter is generated. Concentration simulation samples are generated in a single model process call. One sample; Simultaneously, the job is submitted and executed concurrently through the supercomputing job scheduling system. Each of the aforementioned model invocation processes is independent; Step S2, Covariance matrix dimensionality reduction calculation: Based on the assumption of spatiotemporal separability, the prior error covariance matrix of the state variables is... Decomposed into time covariance matrix With spatial covariance matrix The Kronecker product, i.e. ,in, This represents the Kronecker product operator. The dimension is , The dimension is , For time steps, The number of spatial grid cells; The observation operator matrix Divide into blocks according to the time dimension, represented as follows: Each sub-block The dimension is , The total number of observations; In a shared memory architecture, a multi-threaded parallel computing approach is used, which involves traversing the matrix. For each column, calculate the corresponding matrix product. The column blocks are used directly in subsequent inversion solutions, thus avoiding the need to explicitly store the complete... A matrix, where each thread is responsible for calculation. Zhongyu A corresponding time column block; Step S3, Carbon flux inversion solution: Based on the Bayesian ensemble inversion framework, a cost function containing observation and background terms is constructed, and the observation operator matrix is ​​calculated using the ensemble samples generated in step S1. The linearized expression, and the result obtained using step S2. The matrix is ​​used to solve for the posterior values ​​of the state variables that minimize the cost function, thus completing the high-resolution inversion of urban carbon flux.

2. The method according to claim 1, characterized in that, The serially generated number The value range is greater than 1 and less than or equal to 50.

3. The method according to claim 2, characterized in that, The serially generated number The value of is determined through preliminary experiments and is used to balance the running time and computational resource consumption of a single model process. Its optimized value is 15 to 30.

4. The method according to claim 1, characterized in that, The number of parallel calls The value is dynamically determined based on the total number of available computing nodes and the computing resources of each node, ensuring that the total number of samples in the set is within acceptable limits. It meets the inversion accuracy requirements.

5. The method according to claim 1, characterized in that, The atmospheric transport model is the WRF-Chem model or a derivative thereof, and the modifications to the model code include adding controls. Input parameters and internal loop structure for generating concentration field set samples.

6. The method according to claim 1, characterized in that, The supercomputing job scheduling system is one of Slurm, PBS, or LSF, and the concurrent execution is achieved by combining the batch job script of the selected scheduling system with the MPI parallel environment.

7. The method according to claim 1, characterized in that, The multi-threaded parallel computing is implemented using OpenMP, POSIXThreads, or the multi-threading library built into the programming language.

8. The method according to claim 1, characterized in that, The state variable is urban surface carbon flux, and its dimension is... The observation data was acquired from ground-based and / or satellite sources. Concentration data, of which ground-based observation data are station-based discrete observations with a temporal resolution better than or equal to 1 hour and a spatial resolution better than or equal to kilometers.

9. The method according to claim 1, characterized in that, In step S2, the matrix operation task is assigned to the GPU accelerator for execution by extending the computing interface, so as to further improve computing efficiency.

Images