A method for time series analysis with missing values ​​based on independent components

Abstract

The present invention provides a missing value time series analysis method based on independent components, which belongs to the field of signal processing, and includes the following steps: step S1, obtaining a time series of multi-dimensional observation data, preprocessing it, and then stacking it into an observation data matrix; step S2, decomposing the observation data matrix X, and solving to obtain a whitening component matrix Z; step S3, introducing a separation matrix, further decomposing the whitening component matrix Z, and obtaining mutually independent time components and spatial modes; step S4, using the time components and spatial modes to reconstruct the physical source signal of the time series of the observation data. The present invention is a missing value time series analysis method based on independent components, which is dedicated to the non-interpolation processing of incomplete observation data time series, without pre-interpolation of missing values, and can effectively improve the accuracy of physical signal separation and extraction.

Description

Technical Field

[0001] The invention belongs to the technical field of time series data processing and analysis, and specifically relates to an independent component-based missing value time series analysis method. Background Art

[0002] The observational data sequences in various fields contain a wealth of change signals, such as seasonal signals such as long-term trends, annual cycles, and semi-annual cycles. These signals characterize the behavioral information of different physical processes. However, due to the interference of noise, certain means are usually required to effectively extract the required signals. Common signal analysis methods include principal component analysis, empirical orthogonal functions, singular spectrum analysis, and multi-channel singular spectrum analysis, which can separate signals from noise without relying on any prior knowledge. However, the orthogonal decomposition used by these methods only utilizes the second-order statistical information (variance and covariance) of the observed data, so the components extracted are only uncorrelated rather than independent, causing the reconstructed signal to be largely mixed. Therefore, when multiple physical processes exist at the same time, the components extracted by these methods may still be a mixture of different underlying source signals.

[0003] In theory, the observed time series can be regarded as a linear mixture of multiple independent source signals from different physical processes. Since the source signals and the mixing mode are unknown, signal extraction is actually a blind source separation problem. Independent component analysis (ICA) is an effective method to deal with such problems. It can use high-order statistical information to decompose the mixed signal into multiple statistically independent source signals. Thanks to its powerful blind source separation capability, ICA has been widely used in various research fields, especially in the filtering and signal extraction of GNSS coordinate time series and GRACE time-varying gravity field model.

[0004] However, ICA requires that the observed time series to be processed is complete, and data missing is usually inevitable. Therefore, before performing ICA, missing values ​​must be interpolated to obtain a uniformly sampled observed time series. Commonly used interpolation methods include cubic spline, adjusted maximum likelihood, and iterative interpolation. Although the theoretical perspectives of interpolation have their own merits, no matter which interpolation method is used, certain false information will be introduced, especially when the amount of missing data is large and the fluctuations are frequent, resulting in obvious deviations in the extracted signal. Summary of the invention

[0005] In view of the defects of the prior art, the present invention provides a time series analysis method with missing values ​​based on independent components, which can effectively solve the above problems.

[0006] The technical solution adopted by the present invention is as follows:

[0007] The present invention provides a missing value time series analysis method based on independent components, comprising the following steps:

[0008] Step S1, obtaining a time series of multi-dimensional observation data, preprocessing it, and then stacking it into an observation data matrix;

[0009] Step S2, decomposing the observation data matrix X to obtain the whitening component matrix Z;

[0010] Step S3, introducing a separation matrix to further decompose the whitening component matrix Z to obtain independent time components and spatial modes;

[0011] Step S4, reconstructing the physical source signal of the time series of the observed data using the time component and the spatial mode.

[0012] Preferably, in step S1, the time series of the observed data is expressed as: {x(t i ,j):i=1,2,…,m;j=1,2,…,n}, where n represents the number of observed variables and m represents the number of observed epochs; x(t i ,j) represents the tth i The observation data of the jth observation variable of the observation epoch; the preprocessing includes: non-stationarity test and processing, gross error test and elimination; the observation data matrix is ​​expressed as: X = [x(t i ,1),x(t i ,2),…,x(t i ,n)]; where X is an m×n matrix.

[0013] Preferably, step S2 specifically comprises:

[0014] Step S2.1, determine whether there are missing values ​​in the observation data matrix X, if not, execute step S2.2; if yes, execute step S2.3;

[0015] Step S2.2, using the independent component analysis (ICA) method to obtain a whitening component matrix;

[0016] Step S2.3, using the improved independent component analysis (ICA) method, the whitening component matrix is ​​obtained.

[0017] Preferably, step S2.2 is specifically:

[0018] Step S2.2.1, using formula (1), obtain the covariance matrix B of the observation data matrix X:

[0019] B=X T X (1)

[0020] Among them: the covariance matrix B is an n×n matrix;

[0021] Step S2.2.2, according to the covariance matrix B, use formula (2) to obtain the eigenvector matrix V and the eigenvalue matrix Λ;

[0022] B=VΛV T (2)

[0023] Among them: the eigenvector matrix V is an n×n matrix, and each eigenvector is represented by v k (j), k = 1, 2, ..., n; the eigenvalue matrix Λ is an n × n matrix with eigenvalue λ k , k = 1, 2, ..., n, the eigenvalue matrix Λ is a matrix formed by the eigenvalues ​​of the diagonal elements arranged in descending order;

[0024] Step S2.2.3, using formula (3), solve and obtain the whitening component matrix Z:

[0025] Z=XVΛ -12 (3)

[0026] The whitened component matrix without missing values ​​is obtained by solving this problem.

[0027] Preferably, step S2.3 is specifically:

[0028] Step S2.3.1, the observation data that are not missing in the observation data matrix X are called valid observation data; the formulas of the diagonal elements and non-diagonal elements of the covariance matrix B are established using the valid observation data, thereby obtaining the covariance matrix B;

[0029] Specifically, using formula (4), we can obtain the diagonal elements B(k, k) and off-diagonal elements B(k, j) of the covariance matrix B:

[0030]

[0031] Wherein: the covariance matrix B is an n×n matrix, k=1,2,…,n, j=1,2,…,n, k≠j;

[0032] M k represents the observation data set with no missing k-th observation variable; m k Represents the set M k The number of non-missing observations in ;

[0033] M j represents the set of observation data in which the j-th observation variable is not missing;

[0034] M k ∩M j Indicates M k and M jThe intersection of kj Represents the intersection M k ∩M j The number of non-missing observations;

[0035] Step S2.3.2, based on the obtained covariance matrix B, obtain its eigenvector matrix V and eigenvalue matrix Λ, that is: obtain the eigenvector v k (j) and eigenvalues ​​λk, k = 1, 2, ..., n;

[0036] Step S2.3.3, construct the whitening component expression as follows:

[0037]

[0038] in:

[0039] z k (t i ) represents the tth i The kth whitening component of the observation epoch;

[0040] L i Represents the t i The set of observed variables corresponding to the observation data that are not missing in observation epochs is called the effective observation variable set, L i The number of elements in is 0 to n;

[0041] represents the whitened component of the decomposition of non-missing observation data;

[0042] represents the whitened component of the decomposition of missing observations;

[0043] Step S2.3.4, due to the t i The first observation epoch The observed data x(t i ,j) is missing, therefore, according to the reversibility of the observed data and the whitening component, the expression of the reconstructed missing observed data is obtained:

[0044]

[0045] Step S2.3.5, substitute formula (6) into formula (5) to establish the rank deficiency solution equation for the whitened component:

[0046]

[0047] Equation (7) can be reorganized into a matrix form:

[0048] η(t i )=H(t i )Z(ti ) (8)

[0049] in:

[0050]

[0051] Where: η(t i )、H(t i ) and Z(t i ) are all intermediate matrices; η(t i ) is an intermediate matrix based on non-missing observation data, which can be directly calculated; H(t i ) is an intermediate matrix formed based on eigenvectors and eigenvalues. According to the result of step S2.3.2, H(t i ); therefore, in equation (8), only Z(t i ) is an unknown item to be determined;

[0052] Step S2.3.6, using the minimum norm constraint, solve equation (8) to obtain the complete whitened component matrix.

[0053] Preferably, step S2.3.6 is specifically:

[0054] Step S2.3.6.1, since H(t i ) has a rank equal to the set L i The number of elements in , when there are missing observations, is a rank-deficient situation, and there may be multiple solutions, so the minimum norm condition is introduced to constrain it:

[0055] min:Z T (t i )Z(t i ) (10)

[0056] Step S2.3.6.2, combining equation (8) and constraint (10), construct the objective function Ω as:

[0057] Ω(Z(t i ),α)=Z T (t i )Z(t i )+2α T (H(t i )Z(t i )-η(t i )) (11)

[0058] Where α represents the Lagrange multiplier;

[0059] Step S2.3.6.3, according to the Euler-Lagrange theorem, establish the partial derivative equation:

[0060]

[0061] Step S2.3.6.4, from formula (12) we can get:

[0062] Z(t i )=-H T (t i )α (14)

[0063] Substituting formula (14) into formula (13), we can obtain:

[0064] -H(t i )H T (t i )α=η(t i ) (15)

[0065] Since H(t i ) is the rank-deficient matrix, H(t i )H T (t i ) is also a rank-deficient matrix, so the expression of α is:

[0066] α=-(H(t i )H T (t i )) - η(t i ) (16)

[0067] The superscript “–” indicates pseudo-inverse;

[0068] Substituting equation (16) into equation (14), we get the minimum norm solution of the whitening component matrix:

[0069] Z(t i )=H T (t i )[H(t i )H T (t i )] - η(t i ) (17)

[0070] By solving equation (17), we can obtain the whitening component matrix Z(t i ) is the complete solution.

[0071] Preferably, step S3 specifically comprises:

[0072] Step S3.1. Through Step S2, the whitening component matrix Z obtained by solving is an m×n matrix. The most important and uncorrelated information is mainly concentrated in relatively few components. Therefore, only the first r columns are retained, which can fully represent the components with the greatest variability in the original variables, where r < n, thus obtaining a simplified m×r whitening component matrix, denoted as: whitening component matrix Z0;

[0073] Step S3.2. Since the r column vectors of the whitening component matrix Z0 are only uncorrelated, rather than independent of each other, therefore, a separation matrix W = [w1, w2, …, w r is introduced to further rotate the whitening components to make its column vectors as independent as possible; the decomposition formula of the observed data matrix X is expanded as:

[0074]

[0075] where:

[0076] Λ0 is an r×r eigenvalue matrix, which is the eigenvalue matrix composed of the first r rows and the first r columns of the n×n eigenvalue matrix Λ of the whitening component matrix Z, that is, the eigenvalue matrix of the whitening component matrix Z0;

[0077] is an r×n eigenvector matrix, which is the eigenvector matrix composed of the first r rows of the n×n eigenvector matrix V of the whitening component matrix Z, that is, the eigenvector matrix of the whitening component matrix Z0;

[0078] S is an m×r matrix, denoted as: S = [s1, s2, …, s r = Z0W;

[0079] A T is an r×n matrix, A is an n×r matrix, denoted as:

[0080] where, the column vectors of the matrix S represent independent time components TC; the matrix A is a mixing matrix, and its column vectors represent the spatial modes SP corresponding to the time components. The product of each time component and the corresponding spatial mode is an independent component IC; therefore, by giving a suitable separation matrix W, TC and SP can be uniquely determined, thereby determining the independent component IC; conversely, the separation matrix W can be solved by making the column vectors of the matrix S as independent as possible;

[0081] Step S3.3. Use a fixed-point algorithm based on Newton iteration to solve for the separation matrix W, and then according to the obtained separation matrix W, output the time components and spatial modes.

[0082] Preferably, Step S3.3 is specifically:

[0083] Step S3.3.1, separation matrix W = [w1,w2,…,w r Each column vector in ] is represented as: K , K = 1, 2, ..., r; construct the objective function J G (w K )for:

[0084]

[0085] in:

[0086] E(·) represents the mathematical expectation operator; s Gauss is a Gaussian random variable with zero mean and unit variance;

[0087] G(·) is a non-quadratic function with a continuously differentiable second derivative, and is chosen according to the following conditions:

[0088] If there are both super-Gaussian and sub-Gaussian signals, G(s Gauss )=log2cosh(θ1s Gauss ) / θ1,1≤θ1≤2; θ1 is a constant;

[0089] If only super-Gaussian signals exist, G(s Gauss )=-exp(-θ2s Gauss 2 / 2) / θ2,θ2=1; θ2 is a constant;

[0090] If only sub-Gaussian signals exist, G3(s Gauss )=0.25s Gauss 4 ;

[0091] st represents the constraint condition;

[0092] Step S3.3.2, using the Lagrange multiplier method to solve the optimization problem of equation (19), deriving the iterative solution of the separation matrix W as:

[0093]

[0094] in:

[0095] g(·) is the partial derivative of G(·), l is the number of iterations; g′(·) is the partial derivative of g(·); w K (l+1) and w K (l), represent the values ​​after the l+1th and lth iterations, respectively;

[0096] Step S3.3.3, determine w K Whether the difference after two iterations is less than the preset convergence limit ε:

[0097] ||w K (l+1)-w K (l)||<ε (21)

[0098] If the convergence condition is met, the following formula is used to output a time component s K and spatial mode a K :

[0099]

[0100] This results in a time component and a spatial mode.

[0101] Preferably, step S4 is specifically:

[0102] Step S4.1, reconstruct each independent component IC using the following formula: K :

[0103] IC K =s K a K (twenty three)

[0104] Step S4.2, calculate each independent component IC K Contribution rate CR K :

[0105]

[0106] in:

[0107] ||·|| F represents the Frobenius norm, i.e., F-norm; all independent components are arranged in descending order of contribution rate, and the first d significant independent components are truncated using the F test, while setting the condition set E, s in E K and a K The characteristic conditions of the physical source signal must be met simultaneously;

[0108] Step S4.3, use the following formula to get the physical source signal Signal Recon :

[0109]

[0110] That is, the physical source signal Signal Recon is the sum of the independent components that meet the conditions.

[0111] The present invention provides a method for analyzing time series with missing values ​​based on independent components, which has the following advantages:

[0112] The present invention is a missing value time series analysis method based on independent components, which is dedicated to the non-interpolation processing of incomplete observation data time series. It does not need to pre-interpolate missing values ​​and can effectively improve the accuracy of physical signal separation and extraction. BRIEF DESCRIPTION OF THE DRAWINGS

[0113] Figure 1 A schematic flow chart of a method for analyzing time series with missing values ​​based on independent components provided by the present invention.

[0114] Figure 2 This is a diagram showing the missing ratios of 24 GNSS station sequences in North China in three directions in an embodiment of the present invention;

[0115] Figure 3 A comparison chart of the contribution rates of the first 6 independent components obtained by the two algorithms in the embodiment of the present invention;

[0116] Figure 4 The first MIC (first row) obtained by modified ICA and the first IC (last row) obtained by iterative ICA in the embodiment of the present invention, wherein the top and bottom of each sub-image represent the time component and the spatial mode respectively;

[0117] Figure 5 The second MIC (first row) obtained by modified ICA and the second IC (last row) obtained by iterative ICA in the embodiment of the present invention, wherein the top and bottom of each sub-image represent the time component and the spatial mode respectively;

[0118] Figure 6 This is a comparison diagram of the root mean square error of CMS reconstructed by two algorithms in an embodiment of the present invention. DETAILED DESCRIPTION

[0119] In order to make the technical problems, technical solutions and beneficial effects solved by the present invention more clearly understood, the present invention is further described in detail below in conjunction with the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not used to limit the present invention.

[0120] The purpose of the present invention is to overcome the shortcomings of the background technology and provide a method for extracting signals containing time series with missing values ​​by using the ICA method without any interpolation. The beneficial effects of the present invention compared with the prior art are mainly reflected in that the improved ICA method is based on the principle of "using as much data as possible" and only relies on all available data to complete all steps, giving full play to the characteristics of real data, thereby avoiding the deviation introduced by data interpolation, and will not destroy the structure of the signal, so that the extracted signal has higher accuracy and stronger robustness.

[0121] See also Figure 1 The present invention provides a method for analyzing time series with missing values ​​based on independent components, comprising the following steps:

[0122] Step S1, obtaining a time series of multi-dimensional observation data, preprocessing it, and then stacking it into an observation data matrix;

[0123] Step S2, decomposing the observation data matrix X to obtain the whitening component matrix Z;

[0124] Step S3, introducing a separation matrix to further decompose the whitening component matrix Z to obtain independent time components and spatial modes;

[0125] Step S4, reconstructing the physical source signal of the time series of the observed data using the time component and the spatial mode.

[0126] The present invention is a missing value time series analysis method based on independent components, which is dedicated to the non-interpolation processing of incomplete observation data time series. It does not need to pre-interpolate missing values ​​and can effectively improve the accuracy of physical signal separation and extraction.

[0127] The following is a detailed description of steps S1 to S4:

[0128] Step S1, obtaining a time series of multi-dimensional observation data, preprocessing it, and then stacking it into an observation data matrix;

[0129] Specifically, the time series of observation data is expressed as: {x(t i ,j):i=1,2,…,m;j=1,2,…,n}, where n represents the number of observed variables and m represents the number of observed epochs; x(t i ,j) represents the tth i The observation data of the jth observation variable of the observation epoch;

[0130] In this step, the preprocessing includes: non-stationarity inspection and processing, gross error inspection and elimination, etc.

[0131] The stacked observation data matrix is ​​expressed as: X = [x(t i ,1),x(t i ,2),…,x(t i ,n)]; where X is an m×n matrix, and in general, m≥n.

[0132] Therefore, the observation data matrix X is a matrix formed by the observation data of n observation variables and m observation epochs.

[0133] Step S2, decomposing the observation data matrix X to obtain the whitening component matrix Z;

[0134] Step S2 is specifically as follows:

[0135] Step S2.1, determine whether there are missing values ​​in the observation data matrix X, if not, execute step S2.2; if yes, execute step S2.3;

[0136] Step S2.2, using the independent component analysis (ICA) method to obtain a whitening component matrix;

[0137] Step S2.3, using the improved independent component analysis (ICA) method, the whitening component matrix is ​​obtained.

[0138] It can be seen that in the present invention, different methods for solving the whitening component matrix are used to determine whether the observation data matrix X has missing values.

[0139] Step S2.2, when there are no missing values ​​in the observation data matrix X, the independent component analysis ICA method is used to solve the whitening component matrix. The specific method is:

[0140] Step S2.2.1, since there are no missing values ​​in the observation data matrix X, the covariance matrix B of the observation data matrix X can be directly obtained by using formula (1):

[0141] B=X T X (1)

[0142] Among them: the covariance matrix B is an n×n matrix;

[0143] Step S2.2.2, according to the covariance matrix B, use formula (2) to obtain the eigenvector matrix V and the eigenvalue matrix Λ;

[0144] B=VΛV T (2)

[0145] Among them: the eigenvector matrix V is an n×n matrix, and each eigenvector is represented by v k (j), k = 1, 2, ..., n; the eigenvalue matrix Λ is an n × n matrix with eigenvalue λ k , k = 1, 2, ..., n, the eigenvalue matrix Λ is a matrix formed by the eigenvalues ​​of the diagonal elements arranged in descending order, that is, the diagonal elements are λ1>λ2>...>λ n ;

[0146] Step S2.2.3, using formula (3), solve and obtain the whitening component matrix Z:

[0147] Z=XVΛ -1 / 2 (3)

[0148] The whitened component matrix without missing values ​​is obtained by solving this problem.

[0149] The derivation method of the expression of formula (3) is:

[0150] Since the observation data matrix X can be decomposed into: X = ZΛ 12 V T , thus we get the formula of the whitening component matrix Z in formula (3).

[0151] Step S2.3, using the modified ICA method to obtain the whitening component matrix;

[0152] The main idea of ​​this step is to estimate the covariance matrix using effective observation data, derive the rank deficiency solution equation of the whitened component based on the reversibility of the decomposition process in the time domain and frequency domain, add the minimum norm constraint condition, solve the optimization problem to obtain the complete solution of the whitened component; the specific steps are as follows:

[0153] Step S2.3.1, the observation data that are not missing in the observation data matrix X are called valid observation data; the formulas of the diagonal elements and non-diagonal elements of the covariance matrix B are established using the valid observation data, thereby obtaining the covariance matrix B;

[0154] Specifically, using formula (4), we can obtain the diagonal elements B(k, k) and off-diagonal elements B(k, j) of the covariance matrix B:

[0155]

[0156] Wherein: the covariance matrix B is an n×n matrix, k=1,2,…,n, j=1,2,…,n, k≠j;

[0157] M k represents the observation data set with no missing k-th observation variable; m k Represents the set M k The number of non-missing observations in ;

[0158] M j represents the set of observation data in which the j-th observation variable is not missing;

[0159] M k ∩M j Indicates M k and M j The intersection of kj Represents the intersection M k ∩M j The number of non-missing observations;

[0160] Step S2.3.2, based on the obtained covariance matrix B, obtain its eigenvector matrix V and eigenvalue matrix Λ, that is: obtain the eigenvector v k (j) and the eigenvalue λ k, k=1,2,…,n;

[0161] Step S2.3.3, construct the whitening component expression as follows:

[0162]

[0163] in:

[0164] z k (t i ) represents the tth i The kth whitening component of the observation epoch;

[0165] L i Represents the t i The set of observed variables corresponding to the observation data that are not missing in observation epochs is called the effective set of observed variables. Obviously, L i The number of elements in is 0 to n;

[0166] represents the whitened component of the decomposition of non-missing observation data;

[0167] represents the whitened component of the decomposition of missing observations;

[0168] and The principle is:

[0169] Since the observation data matrix X can be decomposed into: X = ZΛ 12 V T , so the complete form of the matrix form Z of the whitening component written in element form is: Then, the non-missing observation data and the missing observation data are expressed separately, and we get and The expression of , thus obtaining formula (5).

[0170] Step S2.3.4, due to the t i The first observation epoch The observed data x(t i ,j) is missing and cannot be calculated directly. Therefore, according to the reversibility of the observed data and the whitening component, that is: X=ZΛ 12 V T and Z = XVΛ -12 , we get the expression of reconstructed missing observation data:

[0171]

[0172] Step S2.3.5, substitute formula (6) into formula (5) to establish the rank deficiency solution equation for the whitened component:

[0173]

[0174] Equation (7) can be reorganized into a matrix form:

[0175] η(t i )=H(t i )Z(t i ) (8)

[0176] in:

[0177]

[0178] Where: η(t i )、H(t i ) and Z(t i ) are all intermediate matrices; η(t i ) is an intermediate matrix based on non-missing observation data, which can be directly calculated; H(t i ) is an intermediate matrix formed based on eigenvectors and eigenvalues. According to the result of step S2.3.2, H(t i ); therefore, in equation (8), only Z(t i ) is an unknown item to be determined;

[0179] Step S2.3.6, using the minimum norm constraint, solve equation (8) to obtain the complete whitened component matrix.

[0180] Step S2.3.6 is specifically:

[0181] Step S2.3.6.1, since H(t i ) has a rank equal to the set L i The number of elements in , when there are missing observations, is a rank-deficient situation, and there may be multiple solutions, so the minimum norm condition is introduced to constrain it:

[0182] min:Z T (t i )Z(t i )(10)

[0183] Step S2.3.6.2, combining equation (8) and constraint (10), construct the objective function Ω as:

[0184] Ω(Z(t i ),α)=Z T (t i )Z(t i )+2α T (H(t i )Z(t i )-η(ti )) (11)

[0185] Where α represents the Lagrange multiplier;

[0186] Step S2.3.6.3, according to the Euler-Lagrange theorem, establish the partial derivative equation:

[0187]

[0188] Step S2.3.6.4, from formula (12) we can get:

[0189] Z(t i )=-H T (t i )α (14)

[0190] Substituting formula (14) into formula (13), we can obtain:

[0191] -H(t i )H T (t i )α=η(t i ) (15)

[0192] Since H(t i ) is the rank-deficient matrix, H(t i )H T (t i ) is also a rank-deficient matrix, so the expression of α is:

[0193] α=-(H(t i )H T (t i )) - η(t i ) (16)

[0194] The superscript “–” indicates pseudo-inverse;

[0195] Substituting equation (16) into equation (14), we get the minimum norm solution of the whitening component matrix:

[0196] Z(t i )=H T (t i )[H(t i )H T (t i )] - η(t i ) (17)

[0197] By solving equation (17), we can obtain the whitening component matrix Z(t i ) is the complete solution.

[0198] Step S3: Introduce a separation matrix to further decompose the whitened component matrix Z to obtain mutually independent temporal components and spatial modes;

[0199] The main idea of this step is: Based on the non-Gaussianity metric criterion, establish a fixed-point estimation algorithm for the separation matrix, iteratively solve the optimal solution of the separation matrix, realize the further separation of the whitened components, and obtain the temporal components and spatial modes that characterize the observed information.

[0200] Specifically, step S3 is as follows:

[0201] Step S3.1: Through step S2, the obtained whitened component matrix Z is an m×n matrix. The most important and uncorrelated information is mainly concentrated in relatively few components. Therefore, only retaining the first r columns can fully characterize the components with the greatest variability in the original variables, where r < n, thus obtaining a simplified m×r whitened component matrix, denoted as: whitened component matrix Z0;

[0202] Step S3.2: Since the r column vectors of the whitened component matrix Z0 are only uncorrelated, rather than independent of each other, therefore, introduce a separation matrix W = [w1, w2, …, w r to further rotate the whitened components to make its column vectors as independent as possible; at this time, expand the decomposition formula of the observed data matrix X as:

[0203]

[0204] Where:

[0205] Λ0 is an r×r eigenvalue matrix, which is the eigenvalue matrix composed of the first r rows and the first r columns of the n×n eigenvalue matrix Λ of the whitened component matrix Z, that is, the eigenvalue matrix of the whitened component matrix Z0;

[0206] is an r×n eigenvector matrix, which is the eigenvector matrix composed of the first r rows of the n×n eigenvector matrix V of the whitened component matrix Z, that is, the eigenvector matrix of the whitened component matrix Z0;

[0207] S is an m×r matrix, denoted as: S = [s1, s2, …, s r = Z0W;

[0208] A T is an r×n matrix, A is an n×r matrix, denoted as:

[0209] Among them, the column vectors of matrix S represent independent time components TC; matrix A is a mixing matrix, and its column vectors represent the spatial modes SP corresponding to the time components. The product of each time component and the corresponding spatial mode is an independent component IC; therefore, given a suitable separation matrix W, TC and SP can be uniquely determined, thereby determining the independent component IC; conversely, the separation matrix W can be solved by making the column vectors of matrix S as independent as possible; the solution method is shown in step S3.3.

[0210] Step S3.3, considering that the scale of observed data is usually large, the Newton-based iterative fixed point algorithm (FastICA) is used to solve the separation matrix W. The algorithm has a cubic order convergence speed, a computational efficiency that is 10 to 100 times that of the traditional gradient method, and good robustness. Then, according to the solved separation matrix W, the time component and spatial mode are output.

[0211] Step S3.3 is specifically as follows:

[0212] Step S3.3.1, separation matrix W = [w1,w2,…,w r Each column vector in ] is represented as: K , K = 1, 2, ..., r; FastICA measures non-Gaussianity through negative entropy and constructs the objective function J G (w K )for:

[0213]

[0214] in:

[0215] E(·) represents the mathematical expectation operator; s Gauss is a Gaussian random variable with zero mean and unit variance;

[0216] G(·) is a non-quadratic function with a continuously differentiable second derivative, and is chosen according to the following conditions:

[0217] If there are both super-Gaussian and sub-Gaussian signals, G(s Gauss )=log2cosh(θ1s Gauss ) / θ1,1≤θ1≤2; θ1 is a constant;

[0218] If only super-Gaussian signals exist, G(s Gauss )=-exp(-θ2s Gauss 2 / 2) / θ2,θ2=1; θ2 is a constant;

[0219] If only sub-Gaussian signals exist, G3(s Gauss )=0.25s Gauss 4 ;

[0220] st represents the constraint condition;

[0221] Step S3.3.2, using the Lagrange multiplier method to solve the optimization problem of equation (19), deriving the iterative solution of the separation matrix W as:

[0222]

[0223] in:

[0224] g(·) is the partial derivative of G(·), l is the number of iterations; g′(·) is the partial derivative of g(·); w K (l+1) and w K (l), represent the values ​​after the l+1th and lth iterations, respectively;

[0225] Step S3.3.3, determine w K Whether the difference after two iterations is less than the preset convergence limit ε:

[0226] ||w K (l+1)-w K (l)||<ε (21)

[0227] If the convergence condition is met, the following formula is used to output a time component s K and spatial mode a K :

[0228]

[0229] This results in time components and spatial modes.

[0230] Step S4, reconstructing the physical source signal of the time series of the observed data using the time component and the spatial mode.

[0231] Specifically, the independent signal is reconstructed according to the spatiotemporal characteristics and action mechanism of the physical source signal. Step S4 is specifically as follows:

[0232] Step S4.1, reconstruct each independent component IC using the following formula: K :

[0233] IC K =s K a K (twenty three)

[0234] Step S4.2, calculate each independent component IC K Contribution rate CR K :

[0235]

[0236] in:

[0237] ||·|| F represents the Frobenius norm, i.e., F-norm; all independent components are arranged in descending order of contribution rate, and the first d significant independent components are truncated using the F test, while setting the condition set E, s in E K and a K The characteristic conditions of the physical source signal must be met simultaneously;

[0238] Step S4.3, use the following formula to get the physical source signal Signal Recon :

[0239]

[0240] That is, the physical source signal Signal Recon is the sum of the independent components that meet the conditions.

[0241] The present invention provides a method for analyzing time series with missing values ​​based on independent components, which can be summarized as follows:

[0242] For the observation data matrix X with n observation variables and m observation epochs, if X has missing data, the covariance matrix B is calculated using the valid non-missing observation data;

[0243] Taking into account the time domain decomposition (X = ZΛ 1 / 2 V T , Z is the whitening component matrix, Λ and V are the eigenvalue matrix and eigenvector matrix respectively) and frequency domain decomposition (Z = XVΛ -1 / 2 ) process, once the whitened component matrix Z is available, the missing observation data can be reconstructed through the time domain decomposition. Therefore, when the missing data in the frequency domain decomposition is replaced by the data in the time domain decomposition, a rank-deficient system of equations can be obtained to solve the row vectors of the whitened component matrix. In order to solve the rank-deficient equation, the minimum norm criterion is introduced, so that the unique solution of the whitened component matrix can be solved, and this solution is a complete solution.

[0244] Based on the solved whitening component matrix Z, on the one hand, it can be used to reconstruct the complete observation data matrix X, and on the other hand, it can be further used to solve the mixing matrix to achieve the final separation of independent signals.

[0245] An embodiment is described below:

[0246] This embodiment uses the extraction of the common mode signal (CMS) in the GNSS coordinate sequence as an implementation example. CMS is widely considered to be the main source of spatial correlation errors in the GNSS coordinate sequence, and needs to be accurately quantified and separated from the sequence. Previous studies have used interpolation methods to pre-interpolate missing values ​​before using ICA to extract CMS. The embodiment of the present invention uses the modified ICA method of the present invention to compare with the widely used iterative interpolation ICA method (such as Liu B., King M., Dai W., "Common mode error in Antarctic GPS coordinate time-series on its effect on bedrock-uplift estimates" published in 2018). The data used are 24 GNSS reference station coordinate sequences located in North China provided by the China Land State Network, including North (N), East (E), and Up (U) directions, with a time span from January 1, 2011 to November 30, 2019.

[0247] Step S1: Preprocessing of GNSS raw data:

[0248] The non-geophysical offset caused by equipment failure or replacement was corrected according to the station log file; the geophysical offset caused by co-seismic and post-seismic displacement was corrected according to the earthquake record information, and the correction amount was the difference between the average values ​​of the 7 days before and after the offset; the abnormal solutions with a priori errors greater than 100 mm were eliminated, and the gross errors were further detected in combination with the interquartile range statistics, and all eliminated values ​​were marked as missing data; based on the least squares harmonic model, the constant offset term, linear trend term and seasonal term in the sequence were deducted to obtain the residual time series. Figure 2 The data missing ratios of all station time series after preprocessing are shown, with average values ​​of 6.83, 6.82, and 7.13% in the North, East, and Up directions, respectively.

[0249] Step S2: respectively using the modified ICA provided by the present invention and the traditional iterative interpolation iterative ICA to solve the whitening component Z.

[0250] Step S3: Further solve the separation matrix W based on the whitening components solved by the two methods respectively, and obtain each separated time component and the corresponding spatial mode.

[0251] Step S4: Reconstruct the independent components obtained by the two methods and arrange them in descending order according to contribution rate. Figure 3The contribution rates of the first six independent components were compared. The F test showed that the first two independent components were significant at the 95% confidence level. The contribution rates of the first two independent components (MIC) obtained by modified ICA were 17.42% (N), 18.44% (E) and 17.38% (U), while the contribution rates of the first two independent components (IC) obtained by iterative ICA were 16.21% (N), 17.72% (E) and 16.93% (U). Figure 4 and Figure 5 The time components and spatial modes of the first two independent components belonging to the two methods are further given.

[0252] By convention, each spatial mode is normalized by dividing by the maximum absolute value of its elements, and the corresponding time component is scaled by multiplying the normalization factor, so that all spatial modal response values ​​vary from -100% to 100%. In the N, E, and U directions, the average values ​​of the first two spatial modes of MIC are 49% and -31%, -54% and 34%, 45% and 37%, respectively, while the average values ​​of the first two spatial modes of IC are 42% and -30%, -51% and 31%, 56% and 39%, respectively.

[0253] In addition, the spatial modes of the two methods are highly similar. The response direction of the first spatial mode is completely consistent, and the response intensity is distributed in a step-like manner throughout the study area, gradually increasing from south to north, and the larger response values ​​are mainly concentrated in the northeast. Relatively speaking, the intensity of the second spatial mode response is more uniform in the entire area, but its response direction remains completely consistent.

[0254] Step S5: CMS features are defined as follows (referenced in "Spatiotemporal filtering using principal component analysis and Karhunen-Loeve expansion approaches for regional GPS network analysis" by Dong, D., Fang, P., Bock, Y., et al., 2006): (i) the magnitude of the spatial response modes of most stations (>50%) is significant (>25%); (ii) the signs of the spatial response modes of most stations (>90%) are consistent. Figure 4 and Figure 5 By comparison, the first two MICs and ICs meet the CMS characteristics, so they are reconstructed, and then the energy (SP) of the reconstructed CMS is calculated as follows:

[0255]

[0256] Among them, M j and m j They represent the valid data set and epoch number of the jth station, and n is the number of stations. Table 1 gives the energy and operation time comparison of the two algorithms for reconstructing CMS;

[0257] Table 1:

[0258]

[0259] The energies of CMS reconstructed by modified ICA are 0.82mm(N), 1.03mm(E) and 2.77mm(U), which are respectively greater than the energies of CMS reconstructed by iterative ICA, indicating that modified ICA can extract CMS to a greater extent. At the same time, the overall calculation time of modified ICA is shorter, which is mainly attributed to the instability of iterative interpolation, which affects the global efficiency of iterative ICA to a certain extent.

[0260] Step S6: In order to further compare the accuracy of CMS extraction by the two methods under different data missing amounts, the seven stations with the least original missing amounts are selected for simulation experiments.

[0261] First, the CMS extracted from the valid data of 7 stations using two methods was used as the reference signal; then, 5% to 30% (incremental interval was 5%) of the data were randomly deleted, and modified ICA and iterative ICA were performed on the remaining data to extract CMS. The experiment was repeated 200 times for each deletion ratio to avoid the influence of different data deletion positions. The root mean square error (RMSE) was used to evaluate the difference between the reference CMS and the reconstructed CMS under different deletion ratios. The calculation formula is as follows:

[0262]

[0263] The superscripts “M” and “T” represent the experimental results of modified ICA and iterative ICA, respectively. k They represent the reference CMS and the reconstructed CMS after data deletion, n is the number of measuring stations, and K is the number of repeated experiments.

[0264] Figure 6The RMSE values ​​of the two methods for reconstructing CMS are given. As the data deletion ratio increases, the RMSE of the two methods for reconstructing CMS increases accordingly. However, under the same conditions, the RMSE of modified ICA is always smaller than that of iterative ICA, especially in the U direction. When the deletion ratio reaches 30%, the accuracy of modified ICA in extracting CMS is improved by 14.96% (N), 14.75% (E) and 15.67% (U) respectively compared with iterative ICA. Therefore, the method of the present invention can extract the required signal more accurately.

[0265] The above is only a preferred embodiment of the present invention. It should be pointed out that for ordinary technicians in this technical field, several improvements and modifications can be made without departing from the principle of the present invention. These improvements and modifications should also be considered as the scope of protection of the present invention.

Claims

1. A method for time series analysis with missing values ​​based on independent components, characterized in that: It includes the following steps: Step S1, obtain the time series of multi-dimensional observation data, pre-process it, and then stack it into an observation data matrix; the time series of observation data is expressed as: {x(t i ,j):i=1,2,…,m;j=1,2,…,n}, where n represents the number of observed variables and m represents the number of observed epochs; x(t i ,j) represents the tth i The observation data of the jth observation variable of the observation epoch; Step S2: Decompose the observed data matrix X and solve to obtain the whitened component matrix Z; Step S2 is specifically as follows: Step S2.1: Determine whether there are missing values in the observed data matrix X. If not, execute Step S2.2; if so, execute Step S2.3; Step S2.2: Use the independent component analysis (ICA) method to solve and obtain the whitened component matrix; Step S2.3: Use the improved independent component analysis (ICA) method to solve and obtain the whitened component matrix; Step S2.3 is specifically as follows: Step S2.3.1: The observed data in the observed data matrix X that is not missing is called valid observed data; use the valid observed data to establish the formulas for the diagonal elements and non - diagonal elements of the covariance matrix B, thereby obtaining the covariance matrix B; Specifically, use formula (4) to obtain the diagonal element B(k,k) and non - diagonal element B(k,j) of the covariance matrix B: Where: The covariance matrix B is an n×n matrix, k = 1, 2, …, n, j = 1, 2, …, n, k≠j; M k represents the observation data set with no missing k-th observation variable; m k Represents the set M k The number of non-missing observations in ; M j represents the set of observation data in which the j-th observation variable is not missing; M k ∩M j Indicates M k and M j The intersection of kj Represents the intersection M k ∩M j The number of non-missing observations; Step S2.3.2, based on the obtained covariance matrix B, obtain its eigenvector matrix V and eigenvalue matrix Λ, that is: obtain the eigenvector v k (j) and the eigenvalue λ k , k=1,2,…,n; Step S2.3.3: Construct the whitened component expression as follows: Where: z k (t i ) represents the tth i The kth whitening component of the observation epoch; L i Represents the t i The set of observed variables corresponding to the observation data that are not missing in observation epochs is called the effective observation variable set, L i The number of elements in is 0 to n; represents the whitened component of the decomposition of non-missing observation data; represents the whitened component of the decomposition of missing observations; Step S2.3.4, due to the t i The first observation epoch The observed data x(t i ,j) is missing, therefore, according to the reversibility of the observed data and the whitening component, the expression of the reconstructed missing observed data is obtained: Step S2.3.5: Substitute formula (6) into formula (5) to establish the rank - deficiency solution equation for the whitened component: After organizing equation (7), write it in matrix form: η(t i )=H(t i )Z(t i ) (8) Where: Where: η(t i )、H(t i ) and Z(t i ) are all intermediate matrices; η(t i ) is an intermediate matrix based on non-missing observation data, which can be directly calculated; H(t i ) is an intermediate matrix formed based on eigenvectors and eigenvalues. According to the result of step S2.3.2, H(t i ); therefore, in equation (8), only Z(t i ) is an unknown item to be determined; Step S2.3.6: Use the minimum - norm constraint condition to solve equation (8) to obtain the complete whitened component matrix; Step S3: Introduce a separation matrix to further decompose the whitened component matrix Z to obtain mutually independent time components and spatial modes; Step S4: Use the time components and spatial modes to reconstruct the physical source signal of the time series of the observed data.

2. The method for time series analysis with missing values ​​based on independent components according to claim 1, characterized in that: In step S1, the preprocessing includes: non-stationarity test and processing, gross error test and elimination; the observation data matrix is ​​expressed as: X = [x (t i ,1),x(t i ,2),…,x(t i ,n)]; where X is an m×n matrix.

3. The method for time series analysis with missing values ​​based on independent components according to claim 1, characterized in that: Step S2.2 is specifically as follows: Step S2.2.1: Use formula (1) to obtain the covariance matrix B of the observed data matrix X: B=X T X (1) Where: The covariance matrix B is an n×n matrix; Step S2.2.2: According to the covariance matrix B, use formula (2) to obtain the eigenvector matrix V and eigenvalue matrix Λ; B=VΛV T (2) Among them: the eigenvector matrix V is an n×n matrix, and each eigenvector is represented by v k (j), k = 1, 2, ..., n; the eigenvalue matrix Λ is an n × n matrix with eigenvalue λ k , k = 1, 2, ..., n, the eigenvalue matrix Λ is a matrix formed by the eigenvalues ​​of the diagonal elements arranged in descending order; Step S2.2.3: Use formula (3) to solve and obtain the whitened component matrix Z: From the 18th century -1 / 2 (3) Thus, the whitened component matrix in the case of no missing values is solved.

4. The method for analyzing missing values ​​of time series based on independent components according to claim 1, characterized in that: Step S2.3.6 is specifically as follows: Step S2.3.6.1, since H(t i ) has a rank equal to the set L i The number of elements in , when there are missing observations, is a rank-deficient situation, and there may be multiple solutions, so the minimum norm condition is introduced to constrain it: min:Z T (t i )Z(t i ) (10) Step S2.3.6.2: Combine equation (8) and constraint (10) to construct the objective function Ω as: Ω(Z(t i ),α)=Z T (t i )Z(t i )+2α T (H(t i )Z(t i )-η(t i )) (11) Where α represents the Lagrange multiplier; Step S2.3.6.3: According to the Euler - Lagrange theorem, establish the partial - derivative equality: Step S2.3.6.4: From equation (12), it can be obtained that: Z(t i )=-H T (t i )α (14) Substitute equation (14) into equation (13) to get: -H(t i )H T (t i )a=η(t i ) (15) Since H(t i ) is the rank-deficient matrix, H(t i )H T (t i ) is also a rank-deficient matrix, so the expression of α is: α=-(H(t i )H T (t i ))-η(t i ) (16) Where the superscript "–" represents the pseudo - inverse; Substitute equation (16) into equation (14) to obtain the minimum - norm solution of the whitened component matrix as: Z(t i )=H T (t i )[H(t i )H T (t i )] - η(t i ) (17) By solving equation (17), we can obtain the whitening component matrix Z(t i ) is the complete solution.

5. The method for analyzing missing values ​​of time series based on independent components according to claim 2, characterized in that: Step S3 is specifically as follows: Step S3.1: Through Step S2, the obtained whitened component matrix Z is an m×n matrix. The most important and uncorrelated information is mainly concentrated in relatively few components. Therefore, only the first r columns are retained, which can fully represent the components with the largest variability in the original variables, r < n, thereby obtaining a simplified m×r whitened component matrix, denoted as: the whitened component matrix Z0; Step S3.2: Since the r column vectors of the whitening component matrix Z0 are only uncorrelated but not independent of each other, a separation matrix W = [w1, w2, ..., w r ]The whitening component is further rotated to make its column vectors independent of each other; the decomposition of the observation data matrix X is expanded to: Where: Λ0 is an r×r eigenvalue matrix, which is an eigenvalue matrix composed of the first r rows and first r columns of the n×n eigenvalue matrix Λ of the whitened component matrix Z, that is, the eigenvalue matrix of the whitened component matrix Z0; is the r×n eigenvector matrix, which is the eigenvector matrix composed of the first r rows of the n×n eigenvector matrix V of the whitened component matrix Z, that is, the eigenvector matrix of the whitened component matrix Z0; S is an m×r matrix, expressed as: S=[s1,s2,…,s r ] = Z0W; A T is an r×n matrix, and A is an n×r matrix, expressed as: Among them, the column vectors of the matrix S represent the independent time components TC; the matrix A is a mixing matrix, and its column vectors represent the spatial modes SP corresponding to the time components. The product of each time component and the corresponding spatial mode is an independent component IC; therefore, given a suitable separation matrix W, TC and SP can be uniquely determined, thereby determining the independent component IC; conversely, the separation matrix W can be solved by making the column vectors of the matrix S independent; Step S3.3, using the Newton iteration-based fixed point algorithm to solve the separation matrix W, and then outputting the time component and spatial mode according to the solved separation matrix W.

6. The method for analyzing missing values ​​of time series based on independent components according to claim 5, characterized in that: Step S3.3 is specifically as follows: Step S3.3.1, separation matrix W = [w1,w2,…,w r Each column vector in ] is represented as: K , K = 1, 2, ..., r; construct the objective function J G (w K )for: in: E(·) represents the mathematical expectation operator; s Gauss is a Gaussian random variable with zero mean and unit variance; G(·) is a non-quadratic function with a continuously differentiable second derivative, and is chosen according to the following conditions: If there are both super-Gaussian and sub-Gaussian signals, G(s Gauss )=log2cosh(θ1s Gauss ) / θ1,1≤θ1≤2; θ1 is a constant; If only super-Gaussian signals exist, θ2=1; θ2 is a constant; If only sub-Gaussian signals exist, G3(s Gauss )=0.25s Gauss 4 ; st represents the constraint condition; Step S3.3.2, using the Lagrange multiplier method to solve the optimization problem of equation (19), deriving the iterative solution of the separation matrix W as: in: g(·) is the partial derivative of G(·), l is the number of iterations; g′(·) is the partial derivative of g(·); w K (l+1) and w K (l), represent the values ​​after the l+1th and lth iterations, respectively; Step S3.3.3, determine w K Whether the difference after two iterations is less than the preset convergence limit ε: ||w K (l+1)-w K (l)||<ε (21) If the convergence condition is met, the following formula is used to output a time component s K and spatial mode a K : This results in time components and spatial modes.

7. The method for analyzing missing values ​​of time series based on independent components according to claim 6, characterized in that: Step S4 is specifically as follows: Step S4.1, reconstruct each independent component IC using the following formula: K : IC K =s K a K (23) Step S4.2, calculate each independent component IC K Contribution rate CR K : in: ||·|| F represents the Frobenius norm, i.e., F-norm; all independent components are arranged in descending order of contribution rate, and the first d significant independent components are truncated using the F test, while setting the condition set E, s in E K and a K The characteristic conditions of the physical source signal must be met simultaneously; Step S4.3, use the following formula to get the physical source signal Signal Recon : That is, the physical source signal Signal Recon is the sum of the independent components that meet the conditions.

Images