Inertial sensor noise error prediction method

Abstract

The invention discloses an inertial sensor noise error prediction method which comprises the following steps: modeling noise into a transfer function H (z), performing inverse Fourier transform on the transfer function H (z) according to the convolution theorem, calculating to obtain a predicted value y (n), calculating a coefficient a, and substituting the coefficient a into a formula II to obtain the predicted value y (n). According to the method, the error of the inertial sensor is used as a time domain signal for modeling and prediction, and compared with a traditional average value and polynomial interpolation method, the method has the advantages that time sequence signal modeling is added, and the method has good adaptability to interference terms of the sensor.

Description

Technical Field

[0001] The invention relates to a sensor, and in particular to a method for predicting noise errors of an inertial sensor. Background Art

[0002] In the process of using inertial sensors, there are usually a lot of noise and error terms in the measured quantity. The sensor data can be considered as the sum of actual data and noise and error terms. Therefore, there are many estimation algorithms for noise and error terms in the industry, but this goal is difficult to achieve in actual production and R&D. This is because sensor noise and error have many troublesome characteristics. For example, the unknown source of error and noise,

[0003] Temperature, dynamics, micromechanical structure, ADC, etc. will affect the final sensor output, as well as non-repeatability. As the device ages, the accumulation of working hours will change the characteristics of the sensor. Because of these factors, the state of the total error of the sensor always presents a strong nonlinearity. In practical applications, sensor error is not a simple problem. Summary of the invention

[0004] In order to solve the defects of the above-mentioned prior art, the present invention provides an inertial sensor noise error prediction method. The present invention models and predicts the error of the inertial sensor as a time domain signal. Compared with the traditional average value and polynomial interpolation method, the advantage of the present invention is that the modeling of the time series signal is added, and it has better adaptability to the interference terms of the sensor.

[0005] To achieve the above technical objectives, the present invention adopts the following technical solution: a method for predicting noise errors of an inertial sensor, comprising:

[0006] Modeling prediction: Model the noise as a transfer function H(z), and write H(z) as:

[0007]

[0008] According to the convolution theorem, the transfer function H(z) is inversely transformed to calculate the predicted value y(n).

[0009] y(n)=b 0 x(n)+b 1 x(n-1)+…+b N x(nN)+a 1 y(n-1)+a 2 y(n-2)+…+a M y(nM)

[0010] (Formula 2),

[0011] Among them, a and b are filter coefficients, N and M are the lengths of time domain signals;

[0012] In an all-pole filter, N is recorded as 0, and the predicted value y(n) is a function of the coefficient a and the output y;

[0013] Calculate the coefficient a: When the sensor is stationary, obtain a string of digital signals x[n] to form a quasi-diagonal matrix X. Perform conjugate transposition on the diagonal matrix X, and then use the iterative algorithm to calculate the set α of a. The set α is:

[0014]

[0015] Among them, p is the number of algorithm times, k and ε are iterative process variables;

[0016] Substitute the coefficient a into Formula 2 to obtain the predicted value y(n).

[0017] It also includes the dimension increase of the predicted value y(n):

[0018] Obtain a series of time series data X during the sensor’s static phase, expressed as [x(t 1 ) x(t 2 ) … x(t m )], through the mapping relationship S, X'=[x'(t 1 ) x'(t 2 ) … x'(t m )];

[0019] Let data X be represented by matrix X, and X=UΣV H Since H is equivalent to T in real space, X=UΣV T ;

[0020] According to formula 4, formula 5 and formula 6, the matrix XX is calculated. T and the matrix X T X, and the resulting matrix XX T The eigenvalues ​​and eigenvectors of T The eigenvalues ​​and eigenvectors of X, where D is a matrix of dimension m×m or n×n, U is an m×m orthogonal matrix, V is an n×n orthogonal matrix, and Σ is an m×n matrix with values ​​only on the diagonal;

[0021] XX T =UDU T (Formula 4),

[0022] X T X=VDV T (Formula 5)

[0023] D=Σ 2 (Formula 6);

[0024] The matrix XX T The eigenvectors of are merged into a square matrix to obtain the matrix U, and the matrix X T The eigenvectors of X are merged into a square matrix to obtain the matrix V;

[0025] Using matrix U and matrix V, we can calculate X+ according to formula 7.

[0026] PI=VΣ + U(Formula 7),

[0027] Among them, PI is X+;

[0028] According to formula 8, we can get

[0029]

[0030] in, It can improve the prediction dimension and prediction accuracy of the predicted value y(n).

[0031] In summary, the present invention has achieved the following technical effects:

[0032] The advantage of the existing method of taking the average value is that the amount of calculation is very small, and it hardly brings any computational burden to the embedded system. However, the disadvantage is also obvious. A constant cannot predict the time series data at all, and the prediction ability of the unknown error state is weak. The present invention models and predicts the error of the inertial sensor as a time domain signal. Compared with the traditional average value and polynomial interpolation method, the advantage is that the modeling of the time series signal is added, and it has better adaptability to the interference items of the sensor. DETAILED DESCRIPTION

[0033] The present invention is described in further detail below.

[0034] This specific embodiment is merely an explanation of the present invention and is not a limitation of the present invention. After reading this specification, those skilled in the art may make non-creative modifications to the present embodiment as needed. However, as long as they are within the scope of the claims of the present invention, they are protected by the patent law.

[0035] In the description of the present invention, it should be understood that the terms "center", "longitudinal", "lateral", "length", "width", "thickness", "up", "down", "front", "back", "left", "right", "vertical", "horizontal", "top", "bottom", "inside", "outside", "clockwise", "counterclockwise", "axial", "radial", "circumferential" and the like indicating the orientation or position relationship are based on the orientation or position relationship shown, and are only for the convenience of describing the present invention and simplifying the description, and do not indicate or imply that the device or element referred to must have a specific orientation, be constructed and operated in a specific orientation, and therefore should not be understood as limiting the present invention.

[0036] In addition, the terms "first" and "second" are used for descriptive purposes only and should not be understood as indicating or implying relative importance or implicitly indicating the number of the indicated technical features. Therefore, the features defined as "first" and "second" may explicitly or implicitly include one or more of the features. In the description of the present invention, the meaning of "plurality" is two or more, unless otherwise clearly and specifically defined.

[0037] In the present invention, unless otherwise clearly specified and limited, the terms "installed", "connected", "connected", "fixed" and the like should be understood in a broad sense, for example, it can be a fixed connection, a detachable connection, or an integral connection; it can be a mechanical connection or an electrical connection; it can be a direct connection or an indirect connection through an intermediate medium, it can be the internal connection of two elements or the interaction relationship between two elements. For ordinary technicians in this field, the specific meanings of the above terms in the present invention can be understood according to specific circumstances.

[0038] In the present invention, unless otherwise clearly specified and limited, a first feature being "above" or "below" a second feature may mean that the first and second features are in direct contact, or the first and second features are in indirect contact through an intermediate medium. Moreover, a first feature being "above", "above" or "above" a second feature may mean that the first feature is directly above or obliquely above the second feature, or simply means that the first feature is higher in level than the second feature. A first feature being "below", "below" or "below" a second feature may mean that the first feature is directly below or obliquely below the second feature, or simply means that the first feature is lower in level than the second feature.

[0039] Example:

[0040] A method for predicting noise errors of an inertial sensor, comprising:

[0041] Modeling prediction value: Model the noise as a transfer function H(z), and write H(z) as:

[0042]

[0043] According to the convolution theorem, the transfer function H(z) is inversely transformed to calculate the predicted value y(n).

[0044] y(n)=b 0 x(n)+b 1 x(n-1)+…+b N x(nN)+a 1 y(n-1)+a 2 y(n-2)+…+a M y(nM)

[0045] (Formula 2),

[0046] Among them, a and b are filter coefficients, N and M are the lengths of time domain signals;

[0047] In an all-pole filter, N is recorded as 0, and the predicted value y(n) is a function of the coefficient a and the output y;

[0048] Calculate the coefficient a: When the sensor is stationary, obtain a string of digital signals x[n] to form a quasi-diagonal matrix X. Perform conjugate transposition on the diagonal matrix X, and then use the iterative algorithm to calculate the set α of a. The set α is:

[0049]

[0050] Among them, p is the number of algorithm times, k and ε are iterative process variables;

[0051] Substitute the coefficient a into Formula 2 to obtain the predicted value y(n).

[0052] Specifically, in modeling, since low-cost sensors often carry a large number of error terms and noise, these noises are periodic in the time domain, and even different noise sources will produce noises with different characteristics. Therefore, in the frequency domain, the present invention models these noises into a specific transfer function F(s)=Y(s) / X(s), where s is the Laplace operator, s=σ+jω, and after a Laplace transform of the time domain signal, it can be seen that in the s space, the range of ω is [-pi*fs, pi*fs], the range of σ is [-∞, ∞], and the stable domain is σ<0. Therefore, when designing a linear filter, the poles can be placed anywhere in the s domain.

[0053] In a discrete time-invariant system, the frequency ω is uncertain. Under a finite sampling frequency, the oscillator can oscillate at multiples of ω. Therefore, the z transform in the discrete system corresponds to the Laplace transform. In the z transform, the transfer function can be written as H(z)=Y(z) / X(z), where z=ejω.

[0054] Get a series of discrete time domain signals x(n), then the z transform can be written as follows:

[0055]

[0056] The z-change can be viewed as mapping the s-space of the continuous signal into the z-space, with the stable domain within the circle with the origin of the coordinate system.

[0057] The present invention denotes Y(z) and X(z) of the transfer function H(z) as:

[0058]

[0059] but

[0060] The present invention introduces filter coefficients a and b, so Formula 9 is recorded in the form of Formula 1.

[0061] It should be noted that when the numerator and denominator polynomials are zero, there will be zeros and poles of the discrete time invariant system, respectively. By selecting the zero poles, similar discrete systems can be modeled into three models: all pole filter / autoregressive (AR for short), all zero filter / moving average (MA for short), and zero pole filter / autoregressive moving average (ARMA for short). The main function of these filters is to generate time domain data that conforms to frequency domain characteristics. The polynomial coefficients that are finally designed to conform to the transfer function are filter polynomials. The present invention adopts AR model for modeling.

[0062] The present invention predicts the filter in the frequency domain and realizes the filter in the time domain. According to the convolution theorem, the product in the frequency domain can be obtained by Fourier transform of the convolution in the time domain. The same can be obtained by inverse Fourier transform Therefore, the transfer function polynomial and corresponds to the time domain window delayed signal.

[0063] The definition formula of convolution is:

[0064]

[0065] Where h is the filter coefficient, N and M are the lengths of the time domain signals h and x respectively.

[0066] The present invention records the transfer function H(z) in the form of Formula 1. In the frequency domain, the transfer function is inversely Fourier transformed by the convolution theorem and can be written as the sum of two series of convolutions. Thus, the transfer function can be derived into Formula 2:

[0067] y(n)=b 0 x(n)+b 1 x(n-1)+…+b N x(nN)+a 1 y(n-1)+a 2 y(n-2)+…+a M y(nM),

[0068] In two memory spaces of length N and M respectively, the past system input and output data are stored, and then multiplied with the corresponding filter coefficients to obtain the system output at the current moment.

[0069] The present invention also uses linear prediction coefficients (LPC) to perform autoregressive processing on the signal in the signal window to simplify the model of the time series data and predict the subsequent discrete system output.

[0070] In the all-pole filter, in the time domain, a unidirectional channel is obtained to perform convolution calculation on the time domain signal, and the transfer function is recorded as:

[0071]

[0072] in, It is called the synthesis filter. p (z) is called an inverse filter. From its transfer function, it can be seen that the synthetic filter is an all-pole filter with a gain of 1, and its coefficients are the linear prediction coefficients. Because of the all-pole characteristics, the coefficient b in the convolution y(n) will be omitted, leaving only a string of coefficients a and a fixed gain of 1. That is, in formula 2, the present invention only needs to consider coefficient a.

[0073] The process of solving the coefficient a is as follows:

[0074] Get the signal x in the time domain, the fitting residual is r, and a is the equation {a i} opt = argmin[[r(n)] solution. In the algorithm execution phase, in the initialization static phase, a series of fixed-length time-domain sensor signals are sampled, and the single-axis sensor data is used as an example here.

[0075] Get a string of digital signals x[n]. Because of the special conditional judgment, it can be basically considered that the real unbiased sensor output is is 0, so the residual Among the six-axis inertial sensors, the vertical axis accelerometer is an exception. After the acquisition is completed, the obtained data is used to form a quasi-diagonal matrix X, which is written in the following form:

[0076]

[0077] Among them, a is the filter coefficient, and b is the filter coefficient which has been omitted above.

[0078] Solving a is to solve the least squares problem of the linear system Xa=b, through the normal equation X H Xa=X H b, where H represents the conjugate transpose. In this problem, there is no orthogonal signal or quadrature signal, so it is equivalent to the transpose. The linear equations can be written in the form of Yule-Walker equations.

[0079]

[0080] Where R is a Toeplitz matrix, r = [r(1), r(2)...r(p+1) is the autocorrelation prediction coefficient of the time series signal x, This coefficient represents the impact of past moments on subsequent data within the specified interval and in the same time series.

[0081] Finally, using the Levinson-Durbin iterative algorithm, we can calculate the 2 ) is used to calculate the coefficient a. The result of the iterative algorithm is the equation Where p is the number of algorithm iterations, α is the set of coefficients a. k and ε are iterative process variables, and the main loop of the algorithm iteration is:

[0082]

[0083] ε p =-k p+1 β p

[0084]

[0085] In the implementation, according to the sampling frequency of the sensor, the initialization data of fixed length is adopted. After the initialization data acquisition is completed, the above algorithm is used to calculate the lpc and initialize the dsp window to ensure that the window length is consistent with the filter coefficient length, and initialize the FIR filter (FIR class in the program). After the main algorithm is enabled, at each new moment when the data arrives, the dsp window (DigitalSignalWindow class) is slid, the ring buffer (Ringbuffer class) is updated, and it is iterated repeatedly to realize the update of the window and the estimation of the current error term.

[0086] The present invention can combine past information to obtain the predicted value of the latest moment by calculating lpc. However, due to the limitation of the autocorrelation coefficient, the variable y can only be recursively deduced in a single-dimensional space. In order to better predict the sensor error value, the present invention will combine the Koopman theory to upgrade the dimension of the time domain system of the error.

[0087] It also includes the dimension increase of the predicted value y(n):

[0088] Obtain a series of time series data X during the sensor’s static phase, expressed as [x(t 1 ) x(t 2 ) … x(t m )], through the mapping relationship S, X'=[x'(t 1 ) x'(t 2 ) … x'(t m )];

[0089] Let data X be represented by matrix X, and X=UΣV H Since H is equivalent to T in real space, X=UΣV T ;

[0090] According to formula 4, formula 5 and formula 6, the matrix XX is calculated. T and the matrix X T X, and the resulting matrix XX T The eigenvalues ​​and eigenvectors of T The eigenvalues ​​and eigenvectors of X, where D is a matrix of dimension m×m or n×n, U is an m×m orthogonal matrix, V is an n×n orthogonal matrix, and Σ is an m×n matrix with values ​​only on the diagonal;

[0091] XX T =UDU T (Formula 4),

[0092] X T X=VDV T (Formula 5)

[0093] D=Σ 2 (Formula 6);

[0094] The matrix XX T The eigenvectors of are merged into a square matrix to obtain the matrix U, and the matrix X T The eigenvectors of X are merged into a square matrix to obtain the matrix V;

[0095] Using matrix U and matrix V, we can calculate X+ according to formula 7.

[0096] PI=VΣ +U (Formula 7), where PI is X+;

[0097] According to formula 8, we can get in, It can improve the prediction dimension and prediction accuracy of the predicted value y(n),

[0098]

[0099] Specifically, the present invention uses Koopman theory to obtain the function F(x), and establishes a nonlinear mapping relationship y=F(x) between y(n) and F(x).

[0100] According to Koopman theory, if there is a nonlinear equation mapping relationship g is the observation in the state space, S is the mapping to the equation space Γ, and x∈X∈R n If the nonlinear system is finite-dimensional, then there exists an infinite-dimensional linear operator K that satisfies the relation In the equation space F. g∈F(X), because g satisfies the Lipschitz condition, K t :F(X)→F(X) holds, that is, the Koopman operator mapping also exists in the equation space. With these conditions, we can infer that after mapping a state quantity x to the equation space through the mapping relationship g, there exists a finite-dimensional linear operator The observed quantity g is spread over time to approximate the nonlinear motion being addressed.

[0101] In a system, for each eigenvalue λ i , there is an equation Satisfy discrete relations Or in a continuous system Taking out the time partial derivative of the characteristic equation, we can get In the basic properties of linear algebra, there is the equation Av=λv, where v is the eigenvector corresponding to λ.

[0102] Get a string of time series data X in the form of [x(t 1 ) x(t 2 ) … x(t m )], then there is a mapping relation S such that Then X'≈AX is a finite-dimensional estimate of the Koopman operator for the nonlinear system. Assuming that this linear operator can advance the system in time, then in the state vector x∈R m Based on this, there is a linear operator A∈R m×m , in the feature space, there are a finite number of modes that transform the system X in the direction of the eigenvector in the feature space t Convert to Xt+τ If the system is nonlinear, then A is a finite-dimensional estimate of the Koopman operator, giving x k+τ =A Δt x k , that is, A = argmin||X'-AX|| F =X'X + , where F represents the Frobenius norm, + represents the Moore-Penrose pseudo-inverse, and X' is just the time domain data that is overall advanced by a time interval. The solution formula for the Koopman finite estimate is

[0103] According to the principle of singular value decomposition (SVD), the matrix X can be expressed as X=UΣV H , H is equivalent to T in real space, U is an m×m orthogonal matrix, V is an n×n orthogonal matrix, and Σ is an m×n matrix with values ​​only on the diagonal.

[0104] The solution steps are as follows: Because U and V are orthogonal matrices, they satisfy XX T =UDU T , X T X=VDV T , and D = Σ 2 (Write X T =VΣU -1 Substituting the original formula, we can easily prove that in actual operation, the order of solving U and V depends on the dimension of the matrix. In experimental data, the difference between the dimensions of rows and columns is often very large, so the smaller dimension matrix of U and V is calculated first. D is a matrix of m×m or n×n dimensions (depending on whether U or V is solved). According to the matrix diagonalization formula A=PDP -1 , V and U are both orthogonal matrices, so its transpose is equal to the inverse matrix, so U and V are both composed of the corresponding matrix eigenvectors.

[0105] In actual program design, the number of columns of data transmitted to the simulator through the data interface is usually much larger than the number of rows, so U is usually solved first, and the following steps are also derived according to this conclusion.

[0106] First find the matrix XX T =UDU T , and then find the eigenvalue Λ and eigenvector W of the matrix, so that the matrix D can be found. Using the relationship given above, Σ can also be easily calculated. Then merge the obtained eigenvectors into a square matrix to get U. Since the matrix V is too large in this case, we do not completely find the matrix X T X=VDVT In the program, it is pre-set that the matrix X is full rank, and because V is an orthogonal matrix, only m orthogonal vectors are required, and V becomes an n×m matrix.

[0107] Similarly, according to X T X=VDV T Activity matrix V.

[0108] After obtaining the SVD of matrix X, use the formula PI = VΣ + U MP can find the pseudo-inverse, PI is X+, because Σ only has values ​​on the diagonal, so to find the inverse, you only need to take the reciprocal of the value on the diagonal.

[0109] According to the formula *X+ can be calculated

[0110] The present invention utilizes The one-dimensional error prediction is upgraded to multi-dimensional prediction, with higher prediction accuracy.

[0111] In addition, the eigenvalue λ of matrix A is defined as det(A-λI)=0, and the polynomial obtained by finding the determinant is called the characteristic polynomial. The eigenvalues ​​are all the roots of the polynomial. However, it is almost impossible to directly find the roots of the characteristic polynomial in actual implementation, because the roots of high-order polynomials cannot have analytical solutions, and the amount of calculation required to find the determinant increases at a factorial rate, so it is almost impossible even if numerical methods such as Newton's iteration method are used. Therefore, in the program, the matrix must first be converted into Schur Decomposition, S=Q H AQ, Q is a unitary matrix. S is a tridiagonal matrix, which means there are values ​​on the diagonal and upper and lower diagonals, and no values ​​in the lower left part of the matrix. This is also called a Hessenberg matrix. In a Hessenberg matrix, the real eigenvalues ​​appear directly on the main diagonal, while the complex eigenvalues ​​appear as conjugate complex pairs, becoming the roots of the 2x2 matrix on the diagonal. Therefore, if the dimension of the square matrix is ​​odd, there must be a real root.

[0112] The Schur matrix is ​​obtained by using QR decomposition in the program. QR decomposition is defined as A = QR, where Q is a unitary matrix and R is an upper triangular matrix. Using the formula A k+1 =R k Q k =Q k+1 R k+1, iterate a fixed number of times or reach a preset accuracy to get the final Schur form. There are three main methods for QR decomposition, all of which are designed in the program, namely Gram-Schmidt orthogonal method, Givens matrix method, and Householder method.

[0113] Gram-Schmidt splits a matrix into n column vectors, and then orthogonalizes the subsequent vectors with the previous vectors to obtain the final orthogonal matrix. The Geivens rotation matrix obtains the final orthogonal matrix by rotating sub-matrices of different sizes. The Householder method uses the Householder transformation to obtain the final orthogonal matrix. The three methods are similar. The present invention mainly uses the Gram-Schmidt method in the program. At the same time, when the present invention finally calculates the simplified V matrix in the SVD, it does not need to operate on the entire matrix, but only needs to obtain an orthogonal vector with the number of state dimensions.

[0114] (Gram-Schimdt orthogonal vector formula).

[0115] EDMD (Extended Mode Decomposition) is an extension of the traditional DMD algorithm. In traditional DMD, the dimension of the linear operator is fixed. Therefore, when estimating the true value of the sensor, the DMD of the output value of a three-axis sensor can only be three-dimensional, which means that only the influence of three main modes can be integrated. In EDMD, the DMD can be further upgraded to any number of modes, which increases the capacity of the mode library and improves the performance of the predictor. Based on the basic measurement formula of KMD (Koopman Mode Decomposition), where λ is the eigenvalue, is the Koopman characteristic equation, v is the Koopman module, in the Koopman subspace (i.e., the Koopman basis vector characteristic space v), EDMD can infinitely approximate the Koopman operator when the limit of j approaches positive infinity, so compared with DMD, EDMD maps the original state vector to an orthogonal Koopman subspace through the index D

[0116] However, this mapping is not unique, and selecting the mapping relationship D in EDMD is often the most complicated. Common methods include Hermite polynomial, radial kernel function, discrete Fourier modal decomposition, etc. Taking Hermite polynomial as an example, through the mapping relationship D = {ψ 0 ,ψ 1 ,…ψ n}=H i (x)Hj (y), the state vector can be transferred into the Koopman subspace to shift the observations in time. In addition to the traditional kernel function mapping method, the deep neural network embedding algorithm can also effectively estimate the Koopman operator.

[0117] The above description is only a preferred embodiment of the present invention and does not limit the present invention in any form. Any simple modification, equivalent changes and modifications made to the above embodiments based on the technical essence of the present invention are within the scope of the technical solution of the present invention.

Claims

1. A method for predicting noise error of an inertial sensor, characterized in that: include, Modeling prediction: Model the noise as a transfer function H(z), and write H(z) as: According to the convolution theorem, the transfer function H(z) is inversely transformed to calculate the predicted value y(n). y(n)=b0x(n)+b1x(n-1)+…+b N x(n-N)+a1y(n-1)+a2y(n-2)+…+a M y(n-M) (Formula 2), Among them, a and b are filter coefficients, N and M are the lengths of time domain signals; In an all-pole filter, N is recorded as 0, and the predicted value y(n) is a function of the coefficient a and the output y; Calculate the coefficient a: When the sensor is stationary, obtain a string of digital signals x[n] to form a quasi-diagonal matrix X. Perform conjugate transposition on the diagonal matrix X, and then use the iterative algorithm to calculate the set α of a. The set α is: Among them, p is the number of algorithm times, k and ε are iterative process variables; Substitute the coefficient a into Formula 2 to obtain the predicted value y(n).

2. The method for predicting noise error of an inertial sensor according to claim 1, characterized in that: It also includes the dimension increase of the predicted value y(n): Get a series of time series data X during the sensor’s static phase, in the form of [x(t1) x(t2) … x(t m )], through the mapping relationship S, X'=[x'(t1) x'(t2) … x'(t m )]; Let data X be represented by matrix X, and X=UΣV H Since H is equivalent to T in real space, X=UΣV T ; According to formula 4, formula 5 and formula 6, the matrix XX is calculated. T and the matrix X T X, and the resulting matrix XX T The eigenvalues ​​and eigenvectors of T The eigenvalues ​​and eigenvectors of X, where D is a matrix of dimension m×m or n×n, U is an m×m orthogonal matrix, V is an n×n orthogonal matrix, and Σ is an m×n matrix with values ​​only on the diagonal; XX T =UDU T (Formula 4), X T X=VDV T (Formula 5) D=Σ 2 (Formula 6); The matrix XX T The eigenvectors of are merged into a square matrix to obtain the matrix U, and the matrix X T The eigenvectors of X are merged into a square matrix to obtain the matrix V; Using matrix U and matrix V, we can calculate X+ according to formula 7. PI=VΣ + U(Formula 7), Among them, PI is X+; According to formula 8, we can get in, It can improve the prediction dimension and prediction accuracy of the predicted value y(n).

Images