A three-dimensional variational assimilation method for satellite observation data based on channel correlation

Abstract

The application relates to the technical field of data processing, and discloses a satellite observation data three-dimensional variational assimilation method based on channel correlation. The method comprises the following steps: obtaining a channel correlation statistical matrix through statistical analysis and processing; constructing a non-diagonal observation error covariance matrix based on the matrix; decomposing the matrix into covariance sub-matrices by using a matrix block algorithm; quickly solving the observation error weight matrix by using a parallel numerical inversion algorithm; and finally performing variational assimilation calculation and processing on a three-dimensional atmospheric state vector to obtain atmospheric state estimation data that fuses channel correlation. The application solves the technical problems of ignoring channel correlation and low calculation efficiency in the existing satellite observation data three-dimensional variational assimilation method.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of data processing, and in particular to a satellite observation data three-dimensional variational assimilation method based on channel correlation. BACKGROUND

[0002] Satellite observation data three-dimensional variational assimilation is a core technology in a numerical weather prediction system, and the existing technology is mainly based on a traditional diagonalization observation error covariance matrix processing method. This method assumes that the observation errors of different satellite channels are independent of each other, simplifies the observation error covariance matrix into a diagonal matrix form, in which all the non-diagonal elements are set to zero, and only the channel noise variance values on the diagonal are retained. The traditional three-dimensional variational assimilation system uses this simplified processing method to reduce the computational complexity, and performs optimal estimation of the atmospheric state by independently processing the observation data of each channel. This method has been widely used in early numerical weather prediction services.

[0003] However, the diagonalization assumption in the prior art ignores the objective correlation between satellite multi-channel observation data, especially the error correlation relationship between adjacent frequency channels due to the same atmospheric observation path and instrument characteristics. This simplified processing results in a large amount of valuable inter-channel correlation information being discarded, making the assimilation system unable to fully utilize the complementarity and constraint of multi-channel observation data, thereby affecting the accuracy and reliability of the atmospheric state estimation. At the same time, the traditional method faces the problem of low computational efficiency when processing hyperspectral satellite data. When the number of channels increases to hundreds, the storage and inversion of the complete covariance matrix become extremely difficult, and the existing numerical algorithms cannot effectively process large-scale non-diagonal covariance matrices.

[0004] Based on the above analysis, the existing technology also has the problem of fast inversion calculation of the channel correlation covariance matrix, that is, how to effectively reduce the computational complexity of matrix operations while maintaining channel correlation information, and how to design an adaptive block algorithm suitable for the channel correlation structure characteristics to balance the contradiction between computational accuracy and computational efficiency. In addition, the existing technology lacks specialized numerical processing methods for channel correlation characteristics, and cannot perform differential processing according to the correlation strength of different channel groups, resulting in waste of computing resources and numerical stability problems. SUMMARY

[0005] The present application provides a satellite observation data three-dimensional variational assimilation method based on channel correlation, which solves the technical problems of ignoring channel correlation and low computational efficiency in the existing satellite observation data three-dimensional variational assimilation method.

[0006] In a first aspect, the application provides a channel correlation-based three-dimensional variational assimilation method for satellite observation data, which comprises: statistically analyzing and processing channel correlation based on satellite multi-channel observation data to obtain a channel correlation statistical matrix; numerically modeling observation error covariance based on the channel correlation statistical matrix to obtain a non-diagonal observation error covariance matrix; decomposing the non-diagonal observation error covariance matrix by using a matrix block algorithm to obtain a block covariance sub-matrix; quickly inverting the block covariance sub-matrix by using a parallel numerical inversion algorithm to obtain an observation error weight matrix; and performing variational assimilation calculation processing on a three-dimensional atmospheric state vector based on the observation error weight matrix to obtain atmospheric state estimation data fused with channel correlation.

[0007] In a second aspect, the application provides a channel correlation-based three-dimensional variational assimilation system for satellite observation data, which comprises:

[0008] an analysis module configured to statistically analyze and process channel correlation based on satellite multi-channel observation data to obtain a channel correlation statistical matrix;

[0009] a modeling module configured to numerically model observation error covariance based on the channel correlation statistical matrix to obtain a non-diagonal observation error covariance matrix;

[0010] a decomposition module configured to decompose the non-diagonal observation error covariance matrix by using a matrix block algorithm to obtain a block covariance sub-matrix;

[0011] an inversion module configured to quickly invert the block covariance sub-matrix by using a parallel numerical inversion algorithm to obtain an observation error weight matrix;

[0012] a variational module configured to perform variational assimilation calculation processing on a three-dimensional atmospheric state vector based on the observation error weight matrix to obtain atmospheric state estimation data fused with channel correlation.

[0013] In a third aspect, a channel correlation-based three-dimensional variational assimilation device for satellite observation data is provided, which comprises: a memory and at least one processor, the memory storing instructions; and the at least one processor invoking the instructions in the memory to cause the channel correlation-based three-dimensional variational assimilation device for satellite observation data to perform the channel correlation-based three-dimensional variational assimilation method described above.

[0014] In a fourth aspect, a computer readable storage medium is provided, and the computer readable storage medium stores instructions which, when executed on a computer, cause the computer to perform the channel correlation based satellite observation data three-dimensional variational assimilation method described above.

[0015] In the technical scheme provided in the present application, the technical feature of obtaining a channel correlation statistical matrix by statistically analyzing and processing channel correlation of satellite multi-channel observation data effectively solves the problem of ignoring the objective correlation between channels in the traditional method, so that the error correlation relationship between adjacent frequency channels due to sharing the same atmospheric observation path and instrument characteristics is accurately quantified and modeled. The technical feature of obtaining a non-diagonal observation error covariance matrix by numerically modeling the observation error covariance according to the channel correlation statistical matrix breaks the limitations of the traditional diagonalization assumption, re-includes a large amount of discarded inter-channel correlation information into the assimilation framework, and significantly improves the information utilization efficiency of multi-channel observation data. The technical feature of obtaining a block covariance sub-matrix by decomposing the non-diagonal observation error covariance matrix using a matrix block algorithm ingeniously decomposes the large-scale covariance matrix into multiple small-scale sub-matrix blocks, which not only maintains important correlation information but also effectively reduces the computational complexity, laying a foundation for subsequent parallel processing. The technical feature of obtaining an observation error weight matrix by quickly inverting the block covariance sub-matrix through a parallel numerical inversion algorithm breaks through the computational bottleneck faced by the traditional method when dealing with large-scale non-diagonal matrices, and significantly improves the computational efficiency through a block parallel strategy.

[0016] In the field of satellite observation data three-dimensional variational assimilation, the matrix block algorithm adaptively groups according to the strength of channel correlation, so that accurate processing is used in the strong correlation channel group and sparse processing is used in the weak correlation area. This differentiated strategy significantly reduces storage requirements and computational complexity while ensuring numerical accuracy, and is particularly suitable for processing large amounts of channel data generated by modern hyperspectral satellite sensors. The parallel numerical inversion algorithm uses a combination of Cholesky decomposition and sparse matrix inversion, fully utilizes the structural characteristics of the block matrix, and realizes effective parallelization of the computational task, which can obtain significant speedup in a multi-core processor environment. The technical feature of performing variational assimilation calculation on the three-dimensional atmospheric state vector according to the observation error weight matrix enables the final atmospheric state estimation data to fully integrate channel correlation information, and the channel data with strong correlation constrain and supplement each other, which improves the spatial consistency and physical reasonableness of atmospheric temperature and humidity field estimation. This is of great significance for improving the accuracy of numerical weather prediction, especially when dealing with complex temperature and humidity structures in the atmospheric boundary layer and the middle troposphere. BRIEF DESCRIPTION OF DRAWINGS

[0017] In order to more clearly illustrate the technical solutions of the embodiments of the present application, the following will briefly introduce the drawings needed to be used in the embodiment description. Obviously, the drawings in the following description are some embodiments of the present application, and other drawings can be obtained by those skilled in the art without any creative effort based on these drawings.

[0018] Figure 1 An embodiment of the satellite observation data three-dimensional variational assimilation method based on channel correlation in the present application;

[0019] Figure 2 An embodiment of the satellite observation data three-dimensional variational assimilation system based on channel correlation in the present application;

[0020] Figure 3 The structure schematic diagram of the satellite observation data three-dimensional variational assimilation equipment based on channel correlation in the embodiment of the present application. DETAILED DESCRIPTION

[0021] The embodiment of the present application provides a satellite observation data three-dimensional variational assimilation method based on channel correlation. The terms "first", "second", "third", "fourth" and the like (if any) in the specification and claims of the present application and the above-mentioned drawings are used to distinguish similar objects, and do not necessarily have to be used to describe a specific order or sequence. It should be understood that the data thus used can be interchanged under appropriate circumstances, so that the embodiments described herein can be implemented in an order other than that illustrated or described herein. In addition, the terms "include" or "have" and any variations thereof are intended to cover non-exclusive inclusion, for example, a process, method, system, product or device including a series of steps or units does not have to be limited to only those steps or units clearly listed, but can include other steps or units not clearly listed or inherent to these processes, methods, products or devices.

[0022] For the convenience of understanding, the specific process of the embodiment of the present application will be described below. Please refer to Figure 1 An embodiment of the satellite observation data three-dimensional variational assimilation method based on channel correlation in the embodiment of the present application includes:

[0023] Step S101, statistical analysis and processing of channel correlation is performed on satellite multi-channel observation data to obtain a channel correlation statistical matrix;

[0024] Step S102, numerical modeling processing of observation error covariance is performed according to the channel correlation statistical matrix to obtain a non-diagonal observation error covariance matrix;

[0025] Step S103, the non-diagonalized observation error covariance matrix is decomposed by using a matrix block algorithm to obtain a block covariance sub-matrix;

[0026] Step S104, the block covariance sub-matrix is processed by fast inversion through a parallel numerical inversion algorithm to obtain an observation error weight matrix.

[0027] Step S105, the three-dimensional atmospheric state vector is processed by variational assimilation according to the observation error weight matrix to obtain atmospheric state estimation data fused with channel correlation.

[0028] It can be understood that the execution subject of the present application can be a satellite observation data three-dimensional variational assimilation system based on channel correlation, and can also be a terminal or a server, and the specific implementation is not limited herein. The server is taken as an example for description of the embodiments of the present application.

[0029] Specifically, the channel correlation analysis process is started by acquiring satellite radiometer multi-channel brightness temperature observation data. First, the center frequency and half-power bandwidth parameters of each channel are extracted to form a channel spectrum feature parameter set, and then a channel frequency difference matrix is generated based on the frequency interval. The matrix records the frequency distance relationship between channels. The historical observation brightness temperature deviation data is input into a statistical analysis module for sample covariance calculation to generate a channel observation deviation covariance sample matrix. The channel correlation length parameter is obtained by performing channel correlation decay function parameter fitting processing on the channel frequency difference matrix. Finally, the channel correlation statistical matrix is obtained by matrix fusion calculation based on the channel correlation length parameter and the observation deviation covariance sample matrix. The channel correlation length parameter is a key parameter for describing the decay speed of channel correlation with frequency difference, and the channel observation deviation covariance sample matrix reflects the statistical relationship between the errors of the actual observation channels.

[0030] The channel correlation statistical matrix is input into a covariance modeling module for matrix element standardization processing to obtain a standardized correlation matrix. At the same time, the instrument noise equivalent temperature difference parameters of each channel are obtained and a channel noise variance diagonal matrix is constructed. The diagonal elements of the diagonal matrix are the noise variance values of each channel, and the non-diagonal elements are zero. The initial observation error covariance matrix is obtained by matrix multiplication operation based on the standardized correlation matrix and the channel noise variance diagonal matrix. Then, the channel weight function is calculated by vertical integration based on the atmospheric radiation transfer model to obtain a physical correlation correction factor between channels, which considers the influence difference of the atmospheric vertical structure on different channel observations. The physical correlation correction factor between channels is input into the initial observation error covariance matrix for element correction processing to obtain a non-diagonalized observation error covariance matrix. The non-diagonal elements of the matrix contain the correlation information between channels, instead of zero in the traditional method.

[0031] The threshold judgment processing is performed on the element correlation strength of the non-diagonalized observation error covariance matrix. The matrix elements are divided into strong correlation element position index and weak correlation element position index by setting the correlation strength threshold. The matrix row and column are clustered and grouped according to the strong correlation element position index to obtain the channel correlation grouping result. The channels with strong correlation are grouped into the same group. The non-diagonalized observation error covariance matrix is divided into strong correlation submatrix block and weak correlation submatrix block based on the channel correlation grouping result. The elements in the strong correlation submatrix block have strong correlation and need to be kept in the complete structure. The elements in the weak correlation submatrix block have weak correlation and are suitable for sparse processing. The weak correlation submatrix block is input into the sparse algorithm for element zero processing to obtain the sparse weak correlation submatrix. Finally, the matrix reorganization processing is performed on the strong correlation submatrix block and the sparse weak correlation submatrix to obtain the block covariance submatrix.

[0032] The matrix block structure analysis processing is performed on the block covariance submatrix to obtain the dimension information and distribution index of each submatrix block. The dimension information records the number of rows and columns of each submatrix block, and the distribution index records the position of each submatrix block in the original matrix. The Cholesky decomposition processing is performed on the strong correlation submatrix block according to the dimension information and the distribution index to obtain the lower triangular factor matrix. The Cholesky decomposition is a numerical method for decomposing a symmetric positive definite matrix into a lower triangular matrix and its transpose matrix product. The forward substitution and backward substitution operation processing is performed based on the lower triangular factor matrix to obtain the inverse matrix of the strong correlation submatrix block. The forward substitution solves the lower triangular linear equation set, and the backward substitution solves the upper triangular linear equation set. The sparse weak correlation submatrix is input into the sparse matrix inversion algorithm for direct inversion processing to obtain the inverse matrix of the weak correlation submatrix block. The sparse matrix inversion algorithm is specially designed to process the matrix containing a large number of zero elements to reduce the calculation amount. The block reorganization processing is performed on the inverse matrix of the strong correlation submatrix block and the inverse matrix of the weak correlation submatrix block to obtain the observation error weight matrix.

[0033] The background field atmospheric state vector and the multi-channel observation data at the current time are acquired, a three-dimensional variation cost function is initialized and constructed by processing the cost function, the three-dimensional variation cost function includes two main components, i.e., a background field error term and an observation error term. The observation error weight matrix is input into the observation term of the three-dimensional variation cost function to perform weight distribution processing to obtain a weighted observation cost term, the function of the observation error weight matrix is to allocate appropriate weights to the observation data of different channels to reflect the error characteristics and the correlation between channels. The gradient vector of the cost function is obtained by performing gradient calculation processing based on the weighted observation cost term and the background field error covariance term, the gradient vector indicates the direction in which the cost function decreases fastest. The updated atmospheric state vector is obtained by performing conjugate gradient iterative optimization processing on the atmospheric state vector according to the gradient vector, the conjugate gradient method is an effective numerical method for solving large linear equations. The atmospheric state estimation data with channel correlation are obtained by performing convergence test and quality control processing on the updated atmospheric state vector.

[0034] In a specific embodiment, the process of performing step S101 can specifically include the following steps:

[0035] The multi-channel brightness temperature observation data of the satellite radiometer are acquired, the center frequency and the half-power bandwidth parameters of each channel are extracted and processed to obtain a channel spectrum feature parameter set;

[0036] The frequency interval between channels is calculated based on the channel spectrum feature parameter set to obtain a channel frequency difference matrix;

[0037] The historical observation brightness temperature deviation data are input into the statistical analysis module to perform sample covariance calculation processing to obtain a channel observation deviation covariance sample matrix;

[0038] The channel correlation length parameter is obtained by performing parameter fitting processing on the channel correlation decay function according to the channel frequency difference matrix;

[0039] The channel correlation statistical matrix is obtained by performing matrix fusion calculation processing based on the channel correlation length parameter and the observation deviation covariance sample matrix.

[0040] Specifically, the process of acquiring multi-channel brightness temperature observation data of satellite radiometer first receives the raw observation data stream from the satellite data transmission link, which contains the brightness temperature values of each channel and the corresponding instrument calibration information. When extracting the center frequency and half-power bandwidth parameters of each channel, the spectral response function data in the instrument configuration file needs to be read. The center frequency refers to the wavelength or frequency position corresponding to the maximum value of the spectral response function of each channel, and the half-power bandwidth refers to the width of the frequency range corresponding to the value of the spectral response function being half of the maximum value. The channel spectral characteristic parameter set contains the center frequency value, half-power bandwidth value, and shape parameter of the spectral response function of each channel, which describes the sensitivity characteristics of each channel to different frequency electromagnetic radiation.

[0041] When calculating the frequency interval between channels based on the channel spectral characteristic parameter set, the center frequencies of all channels need to be sorted by value, and then the frequency difference between adjacent channels and the frequency distance between any two channels are calculated. The channel frequency difference matrix is a symmetric matrix, where the element in the ith row and jth column represents the absolute value of the center frequency difference between the ith channel and the jth channel. The diagonal elements of this matrix are zero because the frequency difference of each channel with itself is zero. During the calculation process, the correspondence between channel number and frequency size needs to be considered to ensure the accuracy of the frequency difference value.

[0042] When inputting historical observation brightness temperature bias data into the statistical analysis module for sample covariance calculation, the difference between the observation data and the corresponding background field data within a certain time period is first collected as bias samples. The sample covariance calculation process involves de-meaning the bias data, i.e., subtracting the average value of the bias data from each channel, and then calculating the product of the de-meaned bias data between different channels and averaging over time. The element in the ith row and jth column of the channel observation bias covariance sample matrix represents the covariance value between the ith channel and the jth channel observation bias, and this matrix reflects the statistical correlation strength between the observation errors of different channels.

[0043] When fitting the channel correlation decay function according to the channel frequency difference matrix, an exponential decay model is used to describe the variation of channel correlation with frequency difference. The channel correlation decay function usually takes the form of a negative exponential, and the correlation value decreases exponentially with the increase of frequency difference. The parameter fitting process determines the characteristic length parameter in the decay function by least squares method, which controls the speed of correlation decay. The larger the channel correlation length parameter value, the slower the decay of channel correlation, the farther the correlation distance between channels, and vice versa.

[0044] The theoretical correlation model is combined with the actual statistical result when the matrix fusion calculation processing is performed based on the channel correlation length parameter and the observation bias covariance sample matrix. The fusion calculation process first constructs a theoretical correlation matrix according to the channel correlation length parameter and the frequency difference matrix, and then performs weighted average or product operation on the theoretical matrix and the observation bias covariance sample matrix. The inter-channel correlation statistical matrix combines the physical mechanism of frequency correlation and the statistical characteristics of actual observation error, and the matrix element value reflects both the proximity of channel frequencies and the error correlation characteristics in actual observation.

[0045] In a specific embodiment, the process of performing step S102 can specifically include the following steps:

[0046] The inter-channel correlation statistical matrix is input into the covariance modeling module for matrix element standardization processing to obtain a standardized correlation matrix.

[0047] The instrument noise equivalent temperature difference parameters of each channel are obtained, and diagonal matrix construction processing is performed on the noise variance to obtain a channel noise variance diagonal matrix.

[0048] Matrix product operation processing is performed based on the standardized correlation matrix and the channel noise variance diagonal matrix to obtain an initial observation error covariance matrix.

[0049] The channel weight function is vertically integrated according to the atmospheric radiation transmission model to obtain a physical inter-channel correlation correction factor.

[0050] The physical inter-channel correlation correction factor is input into the initial observation error covariance matrix for element correction processing to obtain a non-diagonal observation error covariance matrix.

[0051] Specifically, when the inter-channel correlation statistical matrix is input into the covariance modeling module for matrix element standardization processing, it is first necessary to check whether the diagonal elements of the matrix are 1. If the diagonal elements are not 1, normalization operation is required. The standardization processing is completed by dividing each element of the matrix by the square root product of the corresponding row and column diagonal elements, ensuring that the diagonal elements of the matrix are all 1 and the absolute values of the non-diagonal elements are not more than 1. The covariance modeling module is a calculation unit specially designed for processing multivariate statistical relationships. Its main function is to convert the correlation data into a standardized format for subsequent matrix operations. The standardized correlation matrix maintains the relative relationships of the original correlation structure, while eliminating the influence of dimensional differences on subsequent calculations. The symmetry and positive definiteness of the matrix are maintained in the standardization process.

[0052] The process of obtaining the instrument noise equivalent temperature difference parameters of each channel needs to read the noise characteristic data of each channel from the technical specification document of the satellite instrument. The instrument noise equivalent temperature difference refers to the equivalent brightness temperature change amount corresponding to the instrument noise at a given reference temperature. When constructing the diagonal matrix of the noise variance, the square of the noise equivalent temperature difference value of each channel is filled in as the diagonal element, and all non-diagonal elements are set to zero, indicating that the noise of each channel is independent of each other. The dimension of the channel noise variance diagonal matrix is equal to the number of channels, and the element in the ith row and the ith column of the matrix represents the noise variance value of the ith channel. This matrix describes the error size caused by instrument noise in the observation data of each channel.

[0053] When performing matrix multiplication operation based on the standardized correlation matrix and the channel noise variance diagonal matrix, a triple matrix multiplication calculation method is adopted, i.e., multiplying the square root matrix of the channel noise variance diagonal matrix with the standardized correlation matrix, and then multiplying the square root matrix of the noise variance diagonal matrix. During the matrix multiplication operation process, it is necessary to ensure the matching of the matrix dimensions. The standardized correlation matrix and the noise variance diagonal matrix must have the same number of rows and columns. The initial observation error covariance matrix combines the correlation information between channels and the noise characteristics of each channel. The diagonal elements are equal to the noise variance of each channel, and the non-diagonal elements reflect the covariance relationship of the errors between channels.

[0054] When performing vertical integration calculation on the channel weight function based on the atmospheric radiation transfer model, the spectral response function and the vertical distribution data of atmospheric absorption coefficient of each channel are first obtained. The channel weight function describes the sensitivity of each channel to different atmospheric layers, and its value represents the response strength of the observation of this channel to the temperature and humidity change of a specific atmospheric layer. The vertical integration calculation process multiplies the weight function value of each atmospheric layer by the corresponding layer thickness and sums them up to obtain the comprehensive sensitivity characteristics of each channel to the entire atmospheric column. The physical correlation correction factor between channels is determined by comparing the vertical distribution similarity of different channel weight functions. The more similar the weight function distribution of the channels, the larger the physical correlation correction factor, and vice versa. The atmospheric radiation transfer model is a mathematical model established based on physical principles, which is used to calculate the transmission process of electromagnetic radiation in the atmosphere and the observation characteristics of each channel.

[0055] When the inter-channel physical correlation modification factor is input into the initial observation error covariance matrix for element modification processing, the non-diagonal element of the covariance matrix is adjusted by element-by-element multiplication. The element modification processing process keeps the diagonal elements of the covariance matrix unchanged, and only modifies the physical correlation of the non-diagonal elements. The modified element value is equal to the product of the original covariance value and the corresponding physical correlation modification factor. The non-diagonal observation error covariance matrix contains information of statistical correlation, instrument noise characteristics and physical correlation, and its structure is completely different from the diagonal matrix assumed by the traditional method that the channels are independent.

[0056] In a specific embodiment, the process of performing step S103 can specifically include the following steps:

[0057] The element correlation strength of the non-diagonal observation error covariance matrix is subjected to threshold judgment processing to obtain a strong correlation element position index and a weak correlation element position index.

[0058] The matrix rows and columns are subjected to clustering grouping processing according to the strong correlation element position index to obtain a channel correlation grouping result.

[0059] The non-diagonal observation error covariance matrix is subjected to sub-block division processing based on the channel correlation grouping result to obtain a strong correlation sub-matrix block and a weak correlation sub-matrix block.

[0060] The weak correlation sub-matrix block is input into a sparsification algorithm for element zero processing to obtain a sparsified weak correlation sub-matrix.

[0061] The strong correlation sub-matrix block and the sparsified weak correlation sub-matrix are subjected to matrix reorganization processing to obtain a block covariance sub-matrix.

[0062] Specifically, when the element correlation strength of the non-diagonal observation error covariance matrix is subjected to threshold judgment processing, the absolute values of all non-diagonal elements in the matrix need to be calculated first, and then the correlation strength judgment threshold is determined according to experience or theoretical analysis. The threshold judgment processing process is completed by comparing the size relationship between the absolute value of each non-diagonal element and the preset threshold. When the element absolute value is greater than the threshold, it is classified as a strong correlation element, otherwise it is classified as a weak correlation element. The strong correlation element position index records the row and column coordinate information of all strong correlation elements in the matrix, and the weak correlation element position index records the position information of the weak correlation elements. The selection of the correlation strength threshold directly affects the subsequent matrix block effect. If the threshold is too high, too many elements will be classified as weak correlation, resulting in the loss of important correlation information. If the threshold is too low, the block effect will not be obvious and the computational complexity cannot be reduced.

[0063] In the clustering grouping process of the matrix rows and columns according to the strong correlation element position index, a connectivity analysis method in graph theory is used to identify the channel groups that are correlated with each other. In the clustering grouping process, each channel is regarded as a node in a graph, and the strong correlation is regarded as a connection edge between the nodes. A depth-first search or breadth-first search algorithm is used to find all the connected channel subsets. The channel correlation grouping result contains multiple channel groups. The channels in each group are directly or indirectly strongly correlated with each other, and the channel correlation between different groups is relatively weak. The core idea of the clustering grouping algorithm is to rearrange the original matrix according to the correlation structure, so that the strongly correlated channels are gathered together to form a block structure in the matrix.

[0064] In the sub-block division process of the non-diagonalized observation error covariance matrix based on the channel correlation grouping result, the original matrix is divided into multiple sub-matrix blocks according to the size and position information of each channel group. The sub-block division process needs to maintain the symmetry of the matrix. That is, if the ith channel and the jth channel belong to the same group, the element in the ith row and the jth column of the matrix and the element in the jth row and the ith column of the matrix should be included in the same sub-matrix block. The strongly correlated sub-matrix block corresponds to the covariance relationship between the channels in the same channel group. These sub-matrix blocks maintain the original correlation structure and need to be accurately numerically processed. The weakly correlated sub-matrix block corresponds to the covariance relationship between different channel groups. Since the correlation is weak, these sub-matrix blocks are suitable for sparse processing to reduce the computational burden.

[0065] In the element zeroing process of the weakly correlated sub-matrix block input into the sparse algorithm, the elements to be zeroed are determined according to the preset sparsity requirement or error tolerance. The basic principle of the sparse algorithm is to retain the elements with large absolute values in the matrix and set the elements with small absolute values to zero, thereby reducing the number of non-zero elements while maintaining the main features of the matrix. The element zeroing process needs to consider the symmetry of the matrix. When the element in the ith row and the jth column is zeroed, the element in the jth row and the ith column must also be zeroed to maintain the consistency of the matrix structure. The sparse weakly correlated sub-matrix greatly reduces the number of non-zero elements that need to be stored and calculated while maintaining the basic properties of the original matrix. The choice of sparsity degree needs to balance the computational efficiency and numerical accuracy.

[0066] When performing the matrix reorganization process according to the strongly correlated sub-matrix blocks and the sparsified weakly correlated sub-matrix, each sub-matrix block is reassembled into a complete matrix structure according to the row and column order of the original matrix. The matrix reorganization process needs to ensure that the position of each sub-matrix block is consistent with the position of the corresponding channel in the original matrix, while maintaining the symmetry and positive definiteness of the matrix. The block covariance sub-matrix combines the accurate representation of the strongly correlated region and the sparse representation of the weakly correlated region, which not only maintains important correlation information but also reduces the computational complexity of the matrix. The reorganized matrix has a clear block sparse structure, with the strongly correlated sub-matrix blocks remaining dense and the weakly correlated regions showing sparse distribution.

[0067] In a specific embodiment, the process of inputting the weakly correlated sub-matrix block into the sparsification algorithm for element zeroing processing performed by the execution step can specifically include the following steps:

[0068] Performing size sorting processing on the absolute values of the matrix elements in the weakly correlated sub-matrix block to obtain an element value sorting index;

[0069] Performing sparsity threshold setting processing on the weakly correlated sub-matrix block according to the element value sorting index to obtain a sparsification retention threshold parameter;

[0070] Performing numerical zeroing processing on the elements in the weakly correlated sub-matrix block that are less than the threshold based on the sparsification retention threshold parameter to obtain a preliminary sparsified matrix;

[0071] Inputting the preliminary sparsified matrix into the symmetry maintenance algorithm for matrix symmetry correction processing to obtain a symmetric sparsified matrix;

[0072] Performing positive definiteness inspection and numerical stability adjustment processing on the symmetric sparsified matrix to obtain a sparsified weakly correlated sub-matrix.

[0073] Specifically, when performing size sorting processing on the absolute values of the matrix elements in the weakly correlated sub-matrix block, first, all elements in the sub-matrix block need to be traversed and the absolute values of each element need to be calculated, and then these absolute values are arranged from large to small. The size sorting process is completed using efficient sorting algorithms such as quicksort or mergesort, and the row and column position information of each element in the original matrix needs to be recorded during the sorting process. The element value sorting index contains the numerical value and corresponding matrix position coordinates of each element after sorting, and this index data structure allows the subsequent threshold setting and element zeroing operation to accurately locate the specific matrix element. The purpose of the sorting process is to identify important elements with large numerical values and secondary elements with small numerical values in the matrix, providing a basis for the element retention and discard decision in the sparsification process.

[0074] When performing the sparsity threshold setting process on the weakly correlated submatrix block according to the element value ordering index, a suitable threshold parameter needs to be determined according to the expected sparsification degree and the calculation resource constraint. The sparsity threshold setting process selects a suitable boundary point by analyzing the distribution characteristics of the ordered element values, and the selection of the threshold directly determines the number of non-zero elements finally retained and the sparsity degree of the matrix. The sparsification retention threshold parameter is usually determined based on cumulative importance analysis, that is, elements with greater contribution to the overall characteristics of the matrix are retained, and elements with smaller contribution are discarded. The threshold setting process needs to balance between calculation efficiency and numerical accuracy, and a too high threshold will retain too many elements, resulting in insignificant sparsification effect, and a too low threshold will lose important correlation information and affect assimilation accuracy.

[0075] When performing the numerical zeroing process on the elements smaller than the threshold in the weakly correlated submatrix block based on the sparsification retention threshold parameter, the size relationship between the absolute value of each element in the matrix and the preset threshold is checked one by one. The numerical zeroing process sets all elements with an absolute value smaller than the threshold to zero, while keeping elements with an absolute value greater than or equal to the threshold unchanged, and the processing process needs to maintain the storage structure and index information of the matrix. The preliminary sparsified matrix still maintains the dimension and basic structure of the original matrix after the zeroing process, but most of the elements in it have been set to zero, and the distribution of non-zero elements presents an irregular sparse pattern. The zeroing operation directly affects the storage efficiency of the matrix and the complexity of subsequent calculations, and the storage and operation of the sparse matrix usually use special data structures and algorithms for processing.

[0076] When inputting the preliminary sparsified matrix into the symmetry preserving algorithm for matrix symmetry correction processing, it is necessary to check whether each pair of symmetric elements in the matrix is zero or non-zero at the same time. The core logic of the symmetry preserving algorithm is to ensure that the element in the ith row and jth column of the matrix and the element in the jth row and ith column have the same numerical characteristics, and when an inconsistency is found in the symmetric elements, correction processing is needed. The matrix symmetry correction process usually adopts a relatively conservative strategy, that is, when one of the two symmetric elements is zero, the other is also set to zero, and when both elements are non-zero, their average or smaller value is taken as the corrected value. The symmetric sparsified matrix not only maintains the sparse structure after the correction processing, but also ensures the symmetry of the matrix, which is of great significance for subsequent matrix decomposition and inverse operation.

[0077] In the process of positive definite test and numerical stability adjustment for the symmetric sparse matrix, all eigenvalues of the matrix are calculated to determine whether the matrix is positive definite. The positive definite test is completed by analyzing the signs of the eigenvalues. A matrix is positive definite only if all its eigenvalues are positive. A matrix with zero or negative eigenvalues does not satisfy the positive definite requirement. The numerical stability adjustment includes appropriately amplifying small positive eigenvalues and correcting negative eigenvalues. The adjustment process usually adds small positive numbers to the diagonal elements of the matrix to improve the condition number of the matrix. The sparse weakly related sub-matrix after the positive definite adjustment not only maintains the sparse structure but also satisfies the numerical stability requirement. The matrix can maintain good numerical performance in subsequent matrix inversion and decomposition operations.

[0078] In an embodiment, the process of inputting the preliminary sparse matrix into the matrix symmetry correction algorithm for symmetry modification can include the following steps:

[0079] Transposing the preliminary sparse matrix to obtain a transposed matrix of the preliminary sparse matrix;

[0080] Performing element comparison test based on the preliminary sparse matrix and the transposed matrix to obtain an asymmetric element position index;

[0081] Performing average value calculation on the matrix elements according to the asymmetric element position index to obtain a symmetric correction value;

[0082] Inputting the symmetric correction value into the asymmetric element position of the preliminary sparse matrix for numerical replacement to obtain a modified symmetric matrix;

[0083] Re-inspecting the zero element position and maintaining the sparse structure of the modified symmetric matrix to obtain a symmetric sparse matrix.

[0084] Specifically, when transposing the preliminary sparse matrix, the rows and columns of the matrix need to be interchanged, i.e., the element in the i-th row and the j-th column of the original matrix becomes the element in the j-th row and the i-th column of the transposed matrix. The transposition operation is completed by traversing all elements of the original matrix and redistributing their positions in the new matrix. For a sparse matrix, only non-zero elements need to be processed and their row and column index information needs to be adjusted. The transposed matrix of the preliminary sparse matrix maintains the same sparse structure and number of non-zero elements as the original matrix, but the row and column distribution is mirror changed. The transposed matrix provides a comparison benchmark for subsequent symmetry test. Transposition is a basic operation in linear algebra and is used in satellite observation data processing to test whether the covariance matrix satisfies the symmetry requirement. Because the observation error covariance matrix must be a symmetric matrix in theory to ensure numerical stability.

[0085] When performing element comparison and verification processing based on the preliminary sparse matrix and the transposed matrix, the numerical difference of elements in the corresponding positions of the original matrix and its transposed matrix is compared one by one. The element comparison and verification processing process identifies the asymmetric elements by calculating the difference between the element in the i-th row and j-th column of the original matrix and the element in the i-th row and j-th column of the transposed matrix. When the absolute value of the difference exceeds the preset tolerance range, the element pair is marked as asymmetric. The asymmetric element position index records the specific coordinate information of all detected asymmetric elements in the matrix, including the row number, column number, and corresponding numerical difference size. The comparison and verification algorithm uses a double loop structure to traverse the upper triangular part of the matrix, because the lower triangular part of the symmetric matrix is exactly the same as the upper triangular part, and only half of the elements need to be verified to determine the symmetry of the entire matrix. Special attention needs to be paid to the precision of floating point calculations during the verification process, and relative error rather than absolute error is used to determine whether the elements are symmetric.

[0086] When calculating the average value of the matrix elements according to the asymmetric element position index, the arithmetic mean of each pair of asymmetric elements is calculated as the corrected uniform value. The average value calculation process ensures that the corrected matrix not only maintains the main characteristics of the original data but also meets the symmetry requirement. The choice of average value is a compromise between maintaining data integrity and matrix symmetry. The calculation of the symmetric correction value uses the simple arithmetic mean method, i.e., adding the two values of the asymmetric element pair and dividing by 2. This method can maximize the preservation of original data information while eliminating asymmetry. The calculation process of the correction value also needs to consider the influence of sparsification. When one of the asymmetric element pair is zero, the correction value is usually set to zero to maintain the sparse structure. When both elements are non-zero, their average value is calculated as the correction result.

[0087] When the symmetric correction value is input into the numerical replacement processing of the asymmetric element position in the preliminary sparse matrix, the two elements in the symmetric position of the matrix are updated to have the same value. The numerical replacement processing ensures that the element in the i-th row and j-th column of the matrix has the same value as the element in the j-th row and i-th column. The replacement operation needs to be performed synchronously to avoid data inconsistency. The corrected symmetric matrix maintains the sparse structure while meeting the strict symmetry requirement, and the storage and calculation efficiency of the matrix are considered. Special care needs to be taken in index management during the replacement process to ensure the correctness of the operation. The symmetry correction operation directly affects the numerical stability of subsequent matrix decomposition and inversion operations. Only a strictly symmetric matrix can ensure the correct execution of algorithms such as Cholesky decomposition.

[0088] When rechecking the zero element position and maintaining the sparse structure of the modified symmetric matrix, it is necessary to verify whether the symmetrization modification process introduces new zero elements or changes the original sparse mode. The zero element position rechecking process identifies the positions of zero or near-zero values by traversing all elements of the modified matrix and updates the storage structure and index information of the sparse matrix. The sparse structure maintenance process ensures that the modified matrix still maintains an efficient sparse storage format, avoiding the reduction of storage efficiency caused by symmetrization modification. The number and distribution of non-zero elements need to be recalculated during the process. After the complete modification and verification process of the symmetric sparse matrix, it not only meets the mathematical symmetry requirement but also maintains the sparse advantage in calculation.

[0089] In a specific embodiment, the process of performing step S104 can specifically include the following steps:

[0090] Performing matrix block structure analysis on the block covariance submatrix to obtain dimension information and distribution index of each submatrix block;

[0091] Performing Cholesky decomposition on the strongly correlated submatrix block according to the dimension information and distribution index to obtain a lower triangular factor matrix;

[0092] Performing forward substitution and backward substitution operation based on the lower triangular factor matrix to obtain the inverse matrix of the strongly correlated submatrix block;

[0093] Inputting the sparse weakly correlated submatrix into the sparse matrix inversion algorithm for direct inversion to obtain the inverse matrix of the weakly correlated submatrix block;

[0094] Performing block reorganization according to the inverse matrix of the strongly correlated submatrix block and the inverse matrix of the weakly correlated submatrix block to obtain the observation error weight matrix.

[0095] Specifically, when performing matrix block structure analysis on the block covariance submatrix, it is necessary to scan the entire matrix to identify the distribution of block structure and the boundary position of each submatrix block. The matrix block structure analysis process determines the starting row and column position and the ending row and column position of each submatrix block by checking the aggregation mode of non-zero elements in the matrix. The analysis algorithm uses connectivity detection method to identify dense block area and sparse block area. The dimension information of each submatrix block includes the number of rows and columns of each submatrix block, which is obtained by calculating the difference between the starting position and the ending position of the submatrix block. The distribution index records the exact position coordinates of each submatrix block in the original matrix, including the row and column index of the top-left element and the row and column index of the bottom-right element, which ensures that the subsequent block processing operation can accurately locate each submatrix block. The structure analysis process also needs to verify the independence and integrity of each submatrix block to ensure the rationality of the block scheme and the feasibility of subsequent parallel computation.

[0096] When performing Cholesky decomposition on the strongly correlated sub-matrix block according to the dimension information and the distribution index, it is first verified whether each strongly correlated sub-matrix block satisfies the condition of symmetric positive definiteness, because Cholesky decomposition can only be applied to symmetric positive definite matrices. The Cholesky decomposition process decomposes a symmetric positive definite matrix into the product of a lower triangular matrix and the transpose of the lower triangular matrix, and the decomposition algorithm calculates each element of the lower triangular matrix step by step in rows. The first row and first column element of the lower triangular factor matrix is equal to the square root of the diagonal element of the original matrix, and the subsequent elements are calculated by a recursive relationship. The calculation of each element depends on the values of the elements that have been calculated in the front. Numerical stability checks are required during the decomposition process, and appropriate correction measures need to be taken when non-positive definite matrices or numerical overflow conditions are encountered. The advantage of Cholesky decomposition is its high computational efficiency and good numerical stability, making it particularly suitable for processing covariance matrix inversion problems in satellite observation data. The lower triangular factor matrix obtained by decomposition provides basic data for subsequent substitution operations. When performing forward substitution and backward substitution operations based on the lower triangular factor matrix, the forward substitution algorithm is first used to solve the lower triangular linear equation system, and then the backward substitution algorithm is used to solve the upper triangular linear equation system. The forward substitution operation process starts from the first row of the lower triangular factor matrix and solves row by row. The solution of each unknown is based on the values of the previously calculated unknowns and the corresponding matrix elements. The backward substitution operation process starts from the last row of the transpose of the lower triangular factor matrix and solves in reverse. The calculation order is opposite to that of forward substitution, but the principle is similar. The inverse matrix of the strongly correlated sub-matrix block is obtained by combining the two substitution operations. This inverse matrix maintains the symmetry and positive definiteness of the original matrix. The core idea of substitution is to transform the matrix inversion problem into a linear equation system solving problem, avoiding the numerical instability of direct inversion. The time complexity of the substitution algorithm is relatively low and is suitable for parallel computing implementation.

[0097] When inputting the sparse weakly correlated sub-matrix into the sparse matrix inversion algorithm for direct inversion processing, the special structure of the sparse matrix is used to reduce the computational load and storage requirements. The sparse matrix inversion algorithm is specifically designed for matrices containing a large number of zero elements. The algorithm only performs actual calculations on non-zero elements and skips the processing of zero elements. The direct inversion process uses Gaussian Jordan elimination or LU decomposition, but in the implementation process, sparsity is used to optimize computational efficiency. The algorithm maintains the storage format and index structure of the sparse matrix to avoid unnecessary zero element operations. The inverse matrix of the weakly correlated sub-matrix block may produce fill elements during calculation, i.e., positions that were originally zero become non-zero after inversion. The algorithm needs to dynamically manage these newly generated non-zero elements. The key of the sparse inversion algorithm is to balance the calculation accuracy and efficiency, and to optimize the algorithm performance by estimating the fill pattern and using symbolic decomposition techniques.

[0098] When the inverse matrix of the strongly correlated sub-matrix block and the inverse matrix of the weakly correlated sub-matrix block are subjected to the block reorganization processing, each inverse matrix sub-block is spliced into a complete matrix according to the structure of the original block covariance sub-matrix. The block reorganization processing needs to strictly arrange the positions of each sub-matrix block according to the row and column order and position index of the original matrix, so as to ensure that the structure of the reorganized matrix is consistent with the corresponding relationship of the original matrix. The generation process of the observation error weight matrix involves merging different types of inverse matrix sub-blocks into a unified matrix framework, and the reorganization operation needs to handle the boundary connection and data format unification between different sub-blocks. The reorganized observation error weight matrix inherits the block sparse structure of the original covariance matrix, and the matrix is used to allocate appropriate weights to the observation data of different channels in the three-dimensional variational assimilation calculation.

[0099] In a specific embodiment, the process of performing step S105 can specifically include the following steps:

[0100] Obtaining the background field atmospheric state vector and the multi-channel observation data at the current time, initializing and constructing the cost function to obtain the three-dimensional variational cost function;

[0101] Inputting the observation error weight matrix into the observation term of the three-dimensional variational cost function to perform weight distribution processing to obtain a weighted observation cost term;

[0102] Based on the weighted observation cost term and the background field error covariance term, performing gradient calculation processing to obtain a gradient vector of the cost function;

[0103] According to the gradient vector, performing conjugate gradient iterative optimization processing on the atmospheric state vector to obtain an updated atmospheric state vector;

[0104] Performing convergence inspection and quality control processing on the updated atmospheric state vector to obtain atmospheric state estimation data considering channel correlation.

[0105] Specifically, when acquiring the background field atmospheric state vector and the multi-channel observation data at the current time, first, the prediction result at the previous time from the numerical weather prediction model is read as the background field, which contains the numerical distribution of atmospheric physical quantities such as temperature, humidity, and wind speed on three-dimensional space grid points. The atmospheric state vector arranges all grid point values of the three-dimensional atmospheric field into a one-dimensional vector form in a specific order, and each element of the vector corresponds to the numerical value of a certain physical variable at a specific position and height. The multi-channel observation data at the current time is obtained from a satellite data receiving station, and contains brightness temperature observation values of each channel and corresponding time stamp and geographical position information. When initializing and constructing the cost function, an observation operator needs to be established to connect the atmospheric state vector and the observation data, and the observation operator converts the atmospheric state variables into simulated brightness temperature values of each channel through an atmospheric radiation transfer model. The three-dimensional variational cost function contains two main components, a background field error term and an observation error term, the background field error term measures the deviation between the current state vector and the background field, and the observation error term measures the difference between the simulated observation value and the actual observation value.

[0106]

[0107] wherein represents the numerical value of the three-dimensional variational cost function, which is a scalar function about the atmospheric state vector ; represents the atmospheric state vector to be solved, containing the numerical values of atmospheric variables such as temperature, humidity, and wind speed on all grid points; represents the background field atmospheric state vector, which is derived from the prediction result at the previous time from the numerical prediction model; represents the background field error covariance matrix, which describes the error correlation between different variables and different positions of the background field; represents the inverse matrix of the background field error covariance matrix, which is used as the weight matrix of the background field term; represents the actual observation data vector, which contains satellite brightness temperature observation values of each channel; represents the result of the observation operator acting on the state vector, i.e. the simulated observation value calculated by the atmospheric radiation transfer model; represents the observation error covariance matrix, which corresponds to the observation error weight matrix generated in the foregoing step in the present application; represents the inverse matrix of the observation error covariance matrix, which is used as the weight matrix of the observation term; the superscript T represents the transpose operation of the matrix or vector;

[0108] When the observation error weight matrix is input into the observation term of the three-dimensional variational cost function for weight assignment processing, the observation error weight matrix obtained in the foregoing step is used to replace the traditional diagonalized observation error covariance matrix. The weight assignment processing process uses the observation error weight matrix as the weight coefficient in the observation term, and the non-diagonal elements of the matrix reflect the correlation weight relationship between observation data of different channels. The weighted observation cost term combines the observation-simulation deviation vector and the observation error weight matrix through matrix multiplication operation. Large numerical elements in the weight matrix correspond to high-weight observation data, and small numerical elements correspond to low-weight observation data. The introduction of the observation error weight matrix enables the three-dimensional variational cost function to reasonably utilize the correlation information between channels, and a coordinated and consistent weight assignment is given to a channel combination with strong correlation, while a channel with weak correlation remains relatively independent in weight processing.

[0109] When the gradient calculation process is performed based on the weighted observation cost term and the background field error covariance term, the chain rule is used to calculate the partial derivative of the cost function with respect to the atmospheric state vector. The gradient calculation process includes linearization processing of the observation operator and construction of the adjoint operator. The linearized observation operator describes the influence degree of a small change in the atmospheric state vector on the observation value. The gradient vector of the cost function indicates the direction in which the cost function value decreases most quickly, and each component of the gradient vector corresponds to the partial derivative value of the corresponding element in the state vector. The non-diagonal structure of the observation error weight matrix needs to be specially processed in the gradient calculation. The diagonalization processing in the traditional method is replaced by complete matrix-vector multiplication operation, which increases the calculation complexity but makes more full use of information. The role of the adjoint operator is to transfer the gradient information in the observation space back to the atmospheric state space, ensuring the mathematical correctness and physical reasonableness of the gradient calculation.

[0110] When the conjugate gradient iterative optimization process is performed on the atmospheric state vector according to the gradient vector, the conjugate gradient method is used to solve a large linear equation system to find the minimum point of the cost function. The conjugate gradient iterative optimization process starts from the initial background field state vector, performs line search along the negative gradient direction to determine the optimal step length, then updates the state vector and calculates a new gradient vector. In the iteration process, the concept of conjugate direction is used to accelerate convergence. The new search direction is obtained by linear combination of the current gradient and the previous search direction, and the combination coefficient is determined by the conjugate condition. The updated atmospheric state vector is closer to the minimum point of the cost function in each iteration, and the iteration process continues until the convergence condition is met or the maximum number of iterations is reached. The conjugate gradient algorithm is particularly suitable for processing large-scale sparse linear systems, has low memory requirements and fast convergence speed, and the preconditioning technique of the algorithm further improves the convergence performance.

[0111] When performing convergence test and quality control processing on the updated atmospheric state vector, firstly, it is checked whether the modulus of the gradient vector is less than a preset convergence threshold, and significant reduction of the modulus of the gradient vector indicates that the cost function is close to the minimum value. The convergence test processing also includes checking the variation amplitude of the state vector between consecutive iterations and the downward trend of the cost function value, and simultaneous satisfaction of multiple convergence indicators ensures the reliability of the optimization result. The quality control processing performs physical rationality test on the final atmospheric state vector, including value range check of physical quantities such as temperature and humidity and rationality verification of the vertical gradient. The atmospheric state estimation data fusing the channel correlation comprehensively integrates the background field information and the multi-channel satellite observation information, and enables the observation data of different channels to be complementary and constrained to each other by considering the channel correlation, and the introduction of the channel correlation information in the data fusion process significantly improves the accuracy and consistency of the atmospheric state estimation.

[0112] The above describes the satellite observation data three-dimensional variational assimilation method based on channel correlation in the embodiments of the application, and the following describes the satellite observation data three-dimensional variational assimilation system based on channel correlation in the embodiments of the application. Please refer to Figure 2 One embodiment of the satellite observation data three-dimensional variational assimilation system based on channel correlation in the embodiments of the application includes:

[0113] The analysis module is configured to statistically analyze and process the channel correlation through the satellite multi-channel observation data to obtain a channel correlation statistical matrix.

[0114] The modeling module is configured to numerically model the observation error covariance according to the channel correlation statistical matrix to obtain a non-diagonal observation error covariance matrix.

[0115] The decomposition module is configured to decompose the non-diagonal observation error covariance matrix by using a matrix block algorithm to obtain a block covariance sub-matrix.

[0116] The inversion module is configured to quickly invert the block covariance sub-matrix by using a parallel numerical inversion algorithm to obtain an observation error weight matrix.

[0117] The variational module is configured to perform variational assimilation calculation processing on a three-dimensional atmospheric state vector according to the observation error weight matrix to obtain atmospheric state estimation data fusing the channel correlation.

[0118] The above Figure 2 The satellite observation data three-dimensional variational assimilation system based on channel correlation in the embodiments of the application is described in detail from the perspective of modular functional entities, and the following describes the satellite observation data three-dimensional variational assimilation device based on channel correlation in the embodiments of the application from the perspective of hardware processing.

[0119] Referring toFigure 3 The embodiment of the present application also provides a three-dimensional variational assimilation device based on channel-related satellite observation data, which can be a server, and the internal structure thereof can be as shown in the figure. Figure 3 The three-dimensional variational assimilation device based on channel-related satellite observation data comprises a processor, a memory, a display screen, an input device, a network interface and a database connected through a system bus. The processor of the computer is used to provide computing and control capabilities. The memory of the three-dimensional variational assimilation device based on channel-related satellite observation data comprises a non-volatile storage medium and an internal memory. The non-volatile storage medium stores an operating system, a computer program and a database. The internal memory provides an environment for the operating system and the computer program in the non-volatile storage medium. The database of the three-dimensional variational assimilation device based on channel-related satellite observation data is used to store the corresponding data in the embodiment. The network interface of the three-dimensional variational assimilation device based on channel-related satellite observation data is used to communicate with external terminals through network connection. The computer program is executed by the processor to realize the above method.

[0120] Those skilled in the art can understand that Figure 3 The structure shown in the figure is only a block diagram of part of the structure related to the present application, and does not constitute a limitation on the three-dimensional variational assimilation device based on channel-related satellite observation data to which the present application is applied.

[0121] The present application also provides a computer readable storage medium, which can be a non-volatile computer readable storage medium or a volatile computer readable storage medium, and the computer readable storage medium stores instructions, and when the instructions are run on a computer, the computer executes the steps of the three-dimensional variational assimilation method based on channel-related satellite observation data.

[0122] Those skilled in the art can clearly understand that, for the convenience and brevity of description, the specific working process of the above-mentioned system, system and unit can refer to the corresponding process in the foregoing method embodiments, and will not be repeated here.

[0123] The integrated unit, if implemented in the form of a software function unit and sold or used as an independent product, can be stored in a computer readable storage medium. Based on such understanding, the technical solutions of the present application or the entire or part of the technical solutions that essentially contribute to the prior art can be embodied in the form of a software product. The computer software product is stored in a storage medium and includes a plurality of instructions for causing a channel-related satellite observation data three-dimensional variational assimilation equipment (which can be a personal computer, a server, or a network device, etc.) to execute all or part of the steps of the method described in each embodiment of the present application. The foregoing storage medium includes: a U disk, a mobile hard disk, a read-only memory (ROM), a random access memory (RAM), a magnetic disk or an optical disk, and various program code storage media.

[0124] The above embodiments are only used to illustrate the technical solutions of the present application, but not to limit them; although the present application has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that the technical solutions recorded in the foregoing embodiments can be modified, or some technical features can be replaced by equivalents; and these modifications or replacements do not make the corresponding technical solutions deviate from the spirit and scope of the technical solutions of the embodiments of the present application.

Claims

1. A three-dimensional variational assimilation method for satellite observation data based on channel correlation, characterized in that: The method comprises: Through the satellite multi-channel observation data, the channel correlation is statistically analyzed and processed to obtain the inter-channel correlation statistical matrix; Performing numerical modeling on the observation error covariance according to the inter-channel correlation statistical matrix to obtain a non-diagonalized observation error covariance matrix; The non-diagonalized observation error covariance matrix is ​​decomposed by using a matrix block algorithm to obtain a block covariance submatrix, including: performing threshold judgment processing on the element correlation strength of the non-diagonalized observation error covariance matrix to obtain a strong correlation element position index and a weak correlation element position index; clustering and grouping the matrix rows and columns according to the strong correlation element position index to obtain a channel correlation grouping result; sub-blocking the non-diagonalized observation error covariance matrix based on the channel correlation grouping result to obtain a strong correlation submatrix block and a weak correlation submatrix block; inputting the weak correlation submatrix block into a sparsification algorithm to perform element zeroing processing to obtain a sparse weak correlation submatrix; performing matrix reorganization processing according to the strong correlation submatrix block and the sparse weak correlation submatrix to obtain a block covariance submatrix; The block covariance submatrix is ​​quickly inverted by a parallel numerical inversion algorithm to obtain an observation error weight matrix, including: performing matrix block structure analysis on the block covariance submatrix to obtain dimension information and distribution index of each submatrix block; performing Cholesky decomposition on the strongly correlated submatrix block according to the dimension information and distribution index to obtain a lower triangular factor matrix; performing forward substitution and backward substitution operations based on the lower triangular factor matrix to obtain an inverse matrix of the strongly correlated submatrix block; inputting the sparse weakly correlated submatrix into a sparse matrix inversion algorithm for direct inversion processing to obtain an inverse matrix of the weakly correlated submatrix block; performing block reorganization processing according to the inverse matrix of the strongly correlated submatrix block and the inverse matrix of the weakly correlated submatrix block to obtain an observation error weight matrix; The three-dimensional atmospheric state vector is subjected to variational assimilation calculation processing according to the observation error weight matrix to obtain atmospheric state estimation data fused with channel correlation.

2. The three-dimensional variational assimilation method of satellite observation data based on channel correlation according to claim 1 is characterized in that: The method of performing statistical analysis on channel correlation through satellite multi-channel observation data to obtain an inter-channel correlation statistical matrix includes: Obtain multi-channel brightness temperature observation data from satellite radiometers, extract and process the center frequency and half-power bandwidth parameters of each channel, and obtain a channel spectrum characteristic parameter set; Calculating the frequency intervals between channels based on the channel spectrum characteristic parameter set to obtain a channel frequency difference matrix; The historical observation brightness temperature deviation data is input into the statistical analysis module for sample covariance calculation and processing to obtain the inter-channel observation deviation covariance sample matrix; Performing parameter fitting processing on the channel correlation attenuation function according to the channel frequency difference matrix to obtain a channel correlation length parameter; Matrix fusion calculation processing is performed based on the channel correlation length parameter and the observation deviation covariance sample matrix to obtain an inter-channel correlation statistical matrix.

3. The three-dimensional variational assimilation method of satellite observation data based on channel correlation according to claim 1 is characterized in that: The numerical modeling process is performed on the observation error covariance according to the inter-channel correlation statistical matrix to obtain a non-diagonalized observation error covariance matrix, including: Inputting the inter-channel correlation statistical matrix into the covariance modeling module to perform matrix element standardization processing to obtain a standardized correlation matrix; Obtain the instrument noise equivalent temperature difference parameter of each channel, construct a diagonal matrix of the noise variance, and obtain the channel noise variance diagonal matrix; Performing matrix product operation based on the standardized correlation matrix and the channel noise variance diagonal matrix to obtain an initial observation error covariance matrix; According to the atmospheric radiation transfer model, the channel weight function is vertically integrated and calculated to obtain the correction factor of the physical correlation between channels. The inter-channel physical correlation correction factor is input into the initial observation error covariance matrix for element correction processing to obtain a non-diagonalized observation error covariance matrix.

4. The three-dimensional variational assimilation method of satellite observation data based on channel correlation according to claim 1 is characterized in that: The step of inputting the weakly correlated submatrix block into a sparseness algorithm to perform element zeroing processing to obtain a sparse weakly correlated submatrix comprises: Sorting the absolute values ​​of the matrix elements in the weakly correlated submatrix block to obtain element numerical sorting indexes; Performing sparsity threshold setting processing on the weakly correlated submatrix block according to the element numerical sorting index to obtain a sparsification retention threshold parameter; Based on the sparsification retention threshold parameter, the elements in the weakly correlated submatrix block that are smaller than the threshold are set to zero to obtain a preliminary sparsification matrix; Inputting the preliminary sparsified matrix into a symmetry preserving algorithm to perform matrix symmetry correction processing to obtain a symmetric sparsified matrix; A positive definiteness test and numerical stability adjustment process are performed on the symmetric sparsification matrix to obtain a sparsified weakly correlated submatrix.

5. The three-dimensional variational assimilation method of satellite observation data based on channel correlation according to claim 4 is characterized in that: The step of inputting the preliminary sparsified matrix into a symmetry preserving algorithm to perform matrix symmetry correction processing to obtain a symmetric sparsified matrix comprises: Performing a transpose operation on the preliminary sparsification matrix to obtain a transposed matrix of the preliminary sparsification matrix; Performing element comparison and verification processing based on the preliminary sparsification matrix and the transposed matrix to obtain an asymmetric element position index; Performing average calculation processing on matrix elements according to the asymmetric element position index to obtain a symmetric correction value; Inputting the symmetric correction value into the asymmetric element position of the preliminary sparsification matrix for value replacement processing to obtain a corrected symmetric matrix; The corrected symmetric matrix is ​​subjected to a zero element position recheck and sparse structure preservation processing to obtain a symmetric sparse matrix.

6. The three-dimensional variational assimilation method of satellite observation data based on channel correlation according to claim 1 is characterized in that: The step of performing variational assimilation calculation processing on the three-dimensional atmospheric state vector according to the observation error weight matrix to obtain atmospheric state estimation data fused with channel correlation includes: Obtain the background field atmospheric state vector and the multi-channel observation data at the current moment, initialize and construct the cost function, and obtain a three-dimensional variational cost function; Inputting the observation error weight matrix into the observation term of the three-dimensional variational cost function for weight distribution processing to obtain a weighted observation cost term; Performing gradient calculation based on the weighted observation cost term and the background field error covariance term to obtain a gradient vector of the cost function; performing conjugate gradient iterative optimization processing on the atmospheric state vector according to the gradient vector to obtain an updated atmospheric state vector; The updated atmospheric state vector is subjected to convergence test and quality control processing to obtain atmospheric state estimation data fused with channel correlation.

Images