Parallel kalman filtering systems and methods
By breaking down the Kalman filtering process into parallel steps and using parallel submodules, the problem of low efficiency in traditional Kalman filtering is solved, and efficient computation of multi-state filtering is achieved.
Patent Information
- Application Number
- CN202011358316.X
- Authority / Receiving Office
- CN · China
- Patent Type
- Patents(China)
- Current Assignee / Owner
- Filing Date
- 2020-11-27
- Publication Date
- 2026-01-02
- Estimated Expiration
- 2040-11-27
AI Technical Summary
Traditional Kalman filtering methods are inefficient, cannot perform Kalman filtering for multiple states simultaneously, and consume a lot of computational resources.
The Kalman filtering process is broken down into multiple steps, and the computation is performed in parallel using multiple parallel sub-modules, including modules for transmission matrix generation, observation matrix generation, and state transition calculation. Parallel processing is carried out using ASIC, CPU/GPU cores, or FPGA units.
It improves the operating efficiency of Kalman filtering, enabling Kalman filtering in multiple states simultaneously, thus balancing computational efficiency and resource consumption.
Smart Images

Figure CN113054948B_ABST
Abstract
Description
TECHNICAL FIELD
[0001] The present application relates to the field of communication technology, and particularly relates to a parallel Kalman filtering method for digital signals in the field of communication technology. BACKGROUND
[0002] Filtering theory is a theory and method for estimating the state of a system based on the measurement of observable signals of the system and according to certain filtering criteria. According to the Bayesian theory, the optimal estimation is only an ideal solution, and generally, the analytical form thereof cannot be obtained. In the 1960s, R.E.Kalman proposed a Kalman filtering method, which estimates in the time domain in a recursive manner, and the filtering criterion thereof is the minimum mean square error. When the system is linear and Gaussian, it can obtain the minimum mean square error solution of the recursive Bayesian estimation, and this method overcomes the limitations of the classical Wiener filtering algorithm, is an optimal autoregressive data processing algorithm, and is also the most classical and widely used linear filtering method. This algorithm selects a state space method to design a filter, and can estimate the random process of a multi-dimensional system and a non-stationary system. The Kalman filter has a wide range of applications due to its own advantages, such as recursive operation, adaptability, and foresight, and can optimally estimate a linear dynamic system under random interference. The Kalman filter has been applied to many fields, such as aerospace, communication engineering, artificial intelligence, image processing, rail transportation, and oil exploration. However, the traditional Kalman filtering method has low running efficiency, occupies a large amount of computing resources, and cannot simultaneously implement Kalman filtering of multiple states. SUMMARY
[0003] The technical problem to be solved by the present application is to provide a parallel Kalman filtering method, which can effectively improve the running efficiency of Kalman filtering, can simultaneously implement Kalman filtering of multiple states, and also balances the problem between the computing efficiency and the computing resources, as compared with the traditional Kalman filtering method.
[0004] In order to solve the above technical problem, the present application provides a parallel Kalman filtering method, which includes the following steps: the operation process of Kalman is divided into multiple steps, parallel operation is performed between each step, and each step is simultaneously operated by multiple parallel sub-modules.
[0005] In a preferred embodiment of the present invention, the Kalman operation process is further divided into at least the following steps: S1, generating the transfer matrix A; S2, generating the observation matrix H; S3, calculating the state transition of the system; S4, calculating the observed measurement values; S5, generating the state prediction; S6, generating the prediction covariance P_pre(k); S7, calculating the Kalman gain; S8, calculating the optimal state estimate; S9, updating the state covariance; S10, generating the noise V encountered during the measurement process and the noise W generated during the state transition.
[0006] In a preferred embodiment of the present invention, the invention further includes the following steps:
[0007] (A) Calculate the transmission matrix A: Take A as an n*n dimensional matrix.
[0008] (B) Calculate the observation matrix H: H is an n*n dimensional matrix.
[0009] (C) Calculate the state value of the system: Calculate the state transition of the system using the state transition equation X(k)=A*X(k-1)+L*W(k-1);
[0010] (D) Calculate the observed measurements: using the matrix calculated in step (B). and the k-th state matrix The measurement value Z(K) in the k-th state is obtained by multiplying and adding it to the measurement noise matrix V(K) in the k-th state. The calculation equation for the observation matrix is: Z(k)=H*X(k)+V(k);
[0011] (E) Generate state prediction: Calculate the transfer matrix using step (A). The optimal estimate Xkf(k-1) of the (k-1)th state calculated in step (H) is multiplied by a matrix to obtain the predicted value X_pre(k) of the kth state. Xkf(k-1) is then taken as an n-dimensional column vector. Let X_pre(k) be an n-dimensional column vector. The equation for calculating state prediction is: X_pre(k)=A*Xkf(k-1);
[0012] (F) Generate the predicted covariance P_pre(k): Take P_pre(k) as an n*n dimensional matrix. Its calculation equation is P_pre(k)=A*P(k-1)*A'+Q;
[0013] (G) Calculate the Kalman gain: Take the Kalman gain Kg as an n-dimensional column vector. Its calculation equation is Kg(k)=P_pre(k)*H'*inv(H*P_pre(k)*H'+R);
[0014] (H) Calculate the optimal state estimate: Take the optimal state estimate Its calculation equation is Xkf(k)=X_pre(k)+Kg(k)*(Z(k)-H*X_pre(k));
[0015] (I) Update state covariance: Its calculation equation is P(k)=(I-Kg(k)*H)*P_pre(k);
[0016] (J) and (K) represent the noise encountered during the measurement process and the noise generated during the state transition process, respectively. and The covariances of the two noises are R and R, respectively.
[0017] A parallel Kalman filter system includes a transfer matrix generation module, an observation matrix generation module, a state transition calculation module, an observation measurement calculation module, a state prediction generation module, a prediction covariance generation module, a Kalman gain calculation module, a state optimal estimate calculation module, a state covariance update module, a noise generation module during the measurement process, a noise generation module during the state transition, a process noise covariance generation module, and a measurement noise covariance generation module. The state transition calculation module, observation measurement calculation module, state prediction generation module, prediction covariance generation module, and Kalman gain calculation module are all connected to the transfer matrix generation module and the observation matrix generation module. The state optimal estimate calculation module and the state covariance update module are both connected to the observation matrix generation module. The prediction covariance generation module, Kalman gain calculation module, noise generation module during the measurement process, and noise generation module during the state transition are all connected to the process noise covariance generation module and the measurement noise covariance generation module. The noise generation module during the measurement process and the noise generation module during the state transition are respectively connected to the observation measurement calculation module and the state transition calculation module.
[0018] In a preferred embodiment of the present invention, it further includes (1) a transmission matrix generation module capable of performing the function in step (A);
[0019] (2) The matrix generation module can perform the functions in step (B);
[0020] (3) The state transition calculation module includes multiple independent sub-modules in parallel. Each sub-module that calculates the state transition can complete the full function in step (C) X(k)=A*X(k-1)+L*W(k-1);
[0021] (4) the observation measurement calculation module comprises a plurality of independent sub-modules in parallel, each observation measurement calculation sub-module can complete the complete function Z(k) = H*X(k) + V(k) in step (D);
[0022] (5) the state prediction generation module comprises a plurality of independent sub-modules in parallel, each state prediction generation sub-module can complete the complete function X_pre(k) = A*Xkf(k-1) in step (E);
[0023] (6) the prediction covariance generation module comprises a plurality of independent sub-modules in parallel, each prediction covariance generation sub-module can complete the complete function P_pre(k) = A*P(k-1)*A'+Q in step (F);
[0024] (7) the Kalman gain calculation module comprises a plurality of independent sub-modules in parallel, each Kalman gain calculation sub-module can complete the complete function Kg(k) = P_pre(k)*H'*inv(H*P_pre(k)*H'+R) in step (G);
[0025] (8) the state optimal estimate value calculation module comprises a plurality of independent sub-modules in parallel, each state optimal estimate value calculation sub-module can complete the complete function Xkf(k) = X_pre(k) + Kg(k)*(Z(k)-H*X_pre(k) in step (H);
[0026] (9) the state covariance update module comprises a plurality of independent sub-modules in parallel, each state covariance update sub-module can complete the complete function P(k) = (I-Kg(k)*H)*P_pre(k) in step (I);
[0027] (10) the measurement process noise generation module comprises a plurality of independent sub-modules in parallel, each measurement process noise generation sub-module can complete the complete function in step (J);
[0028] (11) the state transition process noise generation module comprises a plurality of independent sub-modules in parallel, each transmission noise generation sub-module can complete the complete function in step (K).
[0029] In a preferred embodiment of the present application, each module or sub-module is an ASIC functional module, one or more CPU / GPU cores, or one or more FPGA units.
[0030] A parallel Kalman filter sub-module data transmission and operation method comprises the following steps:
[0031] 1) AM calculation generates a transmission matrix Through paths 1_1, 1_2, …, 1_n-1, 1_n to SESM1, SESM2, …, SESMn-1, SESMn, respectively;
[0032] 2) AM calculates the generated transmission matrix Through paths 2_1, 2_2, …, 2_n-1, 2_n to SPESM1, SPESM2, …, SPESMn-1, SPESMn, respectively;
[0033] 3) AM calculates the generated transmission matrix Through paths 3_1, 3_2, …, 3_n-1, 3_n to CPESM1, CPESM2, …, CPESMn-1, CPESMn, respectively;
[0034] 4) HM calculates the generated observation matrix Through paths 4_1, 4_2, …, 4_n-1, 4_n to OESM1, OESM2, …, OESMn-1, OESMn, respectively;
[0035] 5) HM calculates the generated observation matrix Through paths 5_1, 5_2, …, 5_n-1, 5_n to KGESM1, KGESM2, …, KGESMn-1, KGESMn, respectively;
[0036] 6) SESM1 calculates the generated X(1)=system initial state Through path 6_1 to SESM2, SESM2 calculates the generated Through path 6_2 to SESM3, SESMn-1 calculates the generated Through path 6_n-1 to SESMn, SESMn calculates the generated Through path 6_n to SESM1, when N After the end operation;
[0037] 7) WGSM1 calculates the generated Through path 7_1 to SESM2, WGSM2 calculates the generated Through path 7_2 to SESM3, WGSMn-1 calculates the generated Through path 7_n-1 to SESMn, WGSMn calculates the generated Through path 7_n to SESM1, when N After the end operation;
[0038] 8) SEsmi computed by SESM1 SEsm2 computed by SESM2 is transmitted to OESM1 through path 8_1 SEsmn-i computed by SESMn-i is transmitted to OESM2 through path 8_2 SEsmn computed by SESMn is transmitted to OESMn-i through path 8_n-i SEsmn is transmitted to OESMn through path 8_n
[0039] 9) VGsmi computed by VGSM1 VGsm2 computed by VGSM2 is transmitted to OESM1 through path 9_1 VGsmn-i computed by VGSMn-i is transmitted to OESM2 through path 9_2 VGsmn computed by VGSMn is transmitted to OESMn-i through path 9_n-i VGsmn is transmitted to OESMn through path 9_n
[0040] 10) Suesmi computed by SUESM1 Suesm2 computed by SUESM2 is transmitted to SPESM2 through path 10_1 Suesmni computed by SUESMn-i is transmitted to SPESM3 through path 10_2 Suesmn computed by SUESMn is transmitted to SPESMn-i through path 10_n-i Suesmk computed by SUESMk is transmitted to SPESM1 through path 10_n when N < K; when N >= K, the operation is finished
[0041] 11) Cuesmi computed by CUESM1 Cuesm2 computed by CUESM2 is transmitted to CPESM2 through path 11_1 Cuesmni computed by CUESMn-i is transmitted to CPESM3 through path 11_2 Cuesmn computed by CUESMn is transmitted to CPESMn-i through path 11_n-i Cuesmk computed by CUESMk is transmitted to CPESM1 through path 11_n when N < K; when N >= K, the operation is finished
[0042] 12) QGM computed by QGM Through paths 12_1, 12_2, …, 12_n-1, 12_n to CPESM1, CPESM2, …, CPESMn-1, CPESMn respectively;
[0043] 13) CPESM1 calculates generated P_pre(1) = system prediction covariance initial value, which is transmitted to KGESM1 through path 13_1, CPESM2 calculates generated which is transmitted to KGESM2 through path 13_2, CPESMn-1 calculates generated which is transmitted to KGESMn-1 through path 13_n-1, CPESMn calculates generated which is transmitted to KGESMn through path 13_n;
[0044] 14) RGM calculates generated measurement covariance R, which is transmitted to KGESM1, KGESM2, …, KGESMn-1, KGESMn respectively through paths 14_1, 14_2, …, 14_n-1, 14_n;
[0045] 15) SPESM1 calculates generated which is transmitted to SUESM1 through path 15_1, SPESM2 calculates generated which is transmitted to SUESM2 through path 15_2, SPESMn-1 calculates generated which is transmitted to SUESMn-1 through path 15_n-1, SPESMn calculates generated which is transmitted to SUESMn through path 15_n;
[0046] 16) KGESM1 calculates generated Kg which is transmitted to SUESM1 through path 16_1, KGESM2 calculates generated Kg which is transmitted to SUESM2 through path 16_2, KGESMn-1 calculates generated Kg which is transmitted to SUESMn-1 through path 16_n-1, KGESMn calculates generated Kg which is transmitted to SUESMn through path 16_n;
[0047] 17) OESM1 calculates generated which is transmitted to SUESM1 through path 17_1, OESM2 calculates generated which is transmitted to SUESM2 through path 17_2, OESMn-1 calculates generated which is transmitted to SUESMn-1 through path 17_n-1, OESMn calculates generated Transmitted to SUESMn through path 17_n;
[0048] 18) KGESM1 calculates the generated Kg Transmitted to CUESM1 through path 18_1, KGESM2 calculates the generated Kg Transmitted to CUESM2 through path 18_2, KGESMn-1 calculates the generated Kg(n-1) Transmitted to CUESMn-1 through path 18_n-1, KGESMn calculates the generated Kg Transmitted to CUESMn through path 18_n;
[0049] 19) CPESM1 calculates the generated Transmitted to CUESM1 through path 19_1, CPESM2 calculates the generated Transmitted to CUESM2 through path 19_2, CPESMn-1 calculates the generated Transmitted to CUESMn-1 through path 19_n-1, CPESMn calculates the generated Transmitted to CUESMn through path 19_n;
[0050] 20) RGM calculates the generated measurement noise covariance R Transmitted to VGSM1, VGSM2, …, VGSMn-1, VGSMn through paths 20_1, 20_2, …, 20_n-1, 20_n respectively;
[0051] 21) QGM calculates the generated process noise covariance matrix Transmitted to WGSM1, WGSM2, …, WGSMn-1, WGSMn through paths 21_1, 21_2, …, 21_n-1, 21_n respectively.
[0052] In a preferred embodiment of the present application, the number of state points of carrier tracking is 20000, Doppler frequency offset DetaTheta, frequency offset change rate Omega0, and frequency offset change rate change rate Omega1 are selected as the state quantities of the system, the system is a third-order system, detav=1e-3, detaw=1e-3, Ts=1e-3, The empirical formula for calculating the observation matrix by HM is H=[1 Ts / 2 Ts^2 / 6]=[1 5.0e-04 1.67e-07]; and the empirical formula for calculating the observation matrix by AM is The empirical formula for calculating the process noise covariance matrix by QGM is
[0053] The empirical formula for calculating the measurement noise covariance of RGM is R = depth^2 / Ts^4 = 1.0e+06; the measurement noise of VGSM is V = N(0, R); the transmission noise of WGSM is W = N(0, Q); the initial values of the system state variables are set as X(1) = [0 - 20 * e^3 0], the initial values of the system observation variables are Z(1) = [0 00], and the initial values of the system state covariance are... Predicting initial values of covariance The initial value of the optimal state estimate is Xkf (1) = [0 -20*10^(3) 0], and the initial value of the state prediction is X_pre(1) = [0 0 0].
[0054] In a preferred embodiment of the present invention, n is 500, then SESMk completes X(k) = A*X(k-1) + L*W(k-1);
[0055] (A) OESMk completes Z(k)=H*X(k)+V(k);
[0056] (C)CPESMk completes P_pre(k)=A*P(k-1)*A'+Q;
[0057] (D)KGESMk completes Kg(k)=P_pre(k)*H'*inv(H*P_pre(k)*H'+R);
[0058] (E)SUESMk completes Xkf(k)=X_pre(k)+Kg(k)*(Z(k)-H*X_pre(k));
[0059] (F)CUESMk completes P(k)=(I-Kg(k)*H)*P_pre(k).
[0060] In a preferred embodiment of the present invention, the invention further includes the following steps:
[0061] 1) Complete 6 rounds of calculations, generating 3000 system state values;
[0062] 2) Set the initial value of the system state variable Omega1, i.e., X(3, 3000) = 80 * 10^(3), and the initial value of the optimal state estimate Xkf(3, 300) = 80 * 10^(3). Keep other values unchanged, complete one round of calculation, and generate 500 system state values.
[0063] 3) Set the initial value of the system state variables to Omega1, i.e., X(3, 3500) = 0, and the initial value of the optimal state estimate to Xkf(3, 3500) = 0. Keep other values unchanged, complete 11 rounds of calculation, and generate 5500 system state values.
[0064] 4) Set the initial value of the system state quantity Omega1, that is, X(3, 9000) = -80*10^(3), the initial value of the state optimal estimation value Xkf(3, 3500) = -80*10^(3), and other values remain unchanged, complete 1 round of operation, and 500 system state values are generated;
[0065] 5) Set the initial value of the system state quantity Omega1, that is, X(3, 9500) = 0, the initial value of the state optimal estimation value Xkf(3, 3500) = 0, and other values remain unchanged, complete 12 rounds of operation, and 6000 system state values are generated;
[0066] 6) Set the initial value of the system state quantity Omega1, that is, X(3, 15500) = 80*10^(3), the initial value of the state optimal estimation value Xkf(3, 3500) = 80*10^(3), and other values remain unchanged, complete 1 round of operation, and 500 system state values are generated;
[0067] 7) Set the initial value of the system state quantity Omega1, that is, X(3, 16000) = 0, the initial value of the state optimal estimation value Xkf(3, 3500) = 0, and other values remain unchanged, complete 8 rounds of operation, and 4000 system state values are generated.
[0068] In a preferred embodiment of the present application, further comprising.
[0069] The beneficial effects of the present application are:
[0070] The parallel Kalman filtering method of the present application can effectively improve the operation efficiency of Kalman filtering, can simultaneously realize Kalman filtering of multiple states, and balances the problem between calculation efficiency and calculation resources, compared with the traditional Kalman filtering method. BRIEF DESCRIPTION OF DRAWINGS
[0071] Figure 1 It is a structure schematic diagram of the parallel Kalman filtering system in the preferred embodiment of the present application;
[0072] Figure 2 It is a local amplification schematic diagram of the parallel Kalman filtering system in the preferred embodiment of the present application Figure 1 ;
[0073] Figure 3 It is a local amplification schematic diagram of the parallel Kalman filtering system in the preferred embodiment of the present application Figure 2 ;
[0074] Figure 4 It is a schematic diagram of the operation principle of the parallel Kalman filtering system in the preferred embodiment of the present application;
[0075] Figure 5 Fig. 2 is a schematic diagram of data transmission and operation relationship among parallel Kalman filtering sub-modules;
[0076] Figure 6 Fig. 3 is a schematic diagram of data transmission and operation relationship among parallel Kalman filtering sub-modules of 500 sub-modules;
[0077] Figure 7 Fig. 4 is Doppler frequency offset DetaTheta in 20000 system states generated by example calculation;
[0078] Figure 8 Fig. 5 is change rate Omega0 of Doppler frequency offset in 20000 system states generated by example calculation;
[0079] Figure 9 Fig. 6 is change rate Omega1 of change rate of Doppler frequency offset in 20000 system states generated by example calculation;
[0080] Figure 10 Fig. 7 is measurement error and Kalman estimation error of 20000 system states generated by example calculation. DETAILED DESCRIPTION
[0081] The present application will be further described below in conjunction with the drawings and specific embodiments, so that those skilled in the art can better understand the present application and implement it. The embodiments are not intended to limit the present application.
[0082] EMBODIMENT
[0083] The embodiment of the present application discloses a parallel Kalman filtering method, referring to Fig. 1, comprising the following steps: splitting the operation process of Kalman into multiple steps, while making parallel operation between each step, also making simultaneous operation of each step by multiple parallel sub-modules. Compared with the traditional Kalman filtering method, the present application can effectively improve the operation efficiency of Kalman filtering, can simultaneously realize Kalman filtering of multiple states, and also balances the problem between calculation efficiency and calculation resources. Figures 1-4
[0084] When splitting the operation process of Kalman, at least the following steps are included:
[0085] S1, generating transmission matrix A;
[0086] S2, generating observation matrix H;
[0087] S3, calculating state transition of system;
[0088] S4, calculating observation measurement value;
[0089] S5, generating state prediction;
[0090] S6, generating a prediction covariance P_pre(k);
[0091] S7, calculating a Kalman gain;
[0092] S8, calculating a state optimal estimate;
[0093] S9, updating a state covariance;
[0094] S10, generating noise V encountered in a measurement process and noise W generated in a state transition process.
[0095] Specifically, referring to FIG. 1, the method comprises the following steps: Figure 4
[0096] (A) A transmission matrix A generating part is used to calculate a transmission matrix A (here we take A as an n*n dimensional matrix ), the state dimension of the system determines the dimension of A, and the operation process of A and the parameters involved in the operation process are given according to the transmission characteristics of the system itself and the correction of the system in the state transition process according to past experience.
[0097] (B) An observation matrix H generating part is used to calculate an observation matrix H (here we take H as an n*n dimensional matrix ), the state dimension of the system determines the dimension of H, and the operation process of H and the parameters involved in the operation process are given according to the to-be-observed quantity of the measurement process (or observation process) and the correction of the measurement system according to past experience.
[0098] (C) A state transition calculation part is used to calculate the state value of the system. The calculation process of the state transition of the system is carried out according to the state transition equation X(k) = A*X(k-1) + L*W(k-1), wherein X (here we take X as an n-dimensional column vector ) is the state matrix of the system, k indicates the kth state X(k) (here we take X(k) as an n-dimensional column vector ) in the state transition process, k-1 indicates the (k-1)th state X(k-1) (here we take X(k-1) as an n-dimensional column vector ) in the state transition process. L (here we take L as an n*n dimensional matrix ) is a parameter matrix for adjusting the noise in the state transition process, the dimension of which is determined by the state dimension of the system, and W (here we take W as an n-dimensional column vector ) will be described in part (K). The state transition equation X(k) = A*X(k-1) + L*W(k-1) can be expressed in matrix form as,
[0099] That is
[0100]
[0101]
[0102] (D) The calculation of the measurement value is the multiplication of the matrix calculated in (B) and the kth state matrix and the measurement noise matrix V(k) in the kth state (here we take V(k) as an n-dimensional column vector ) to get the measurement value Z(k) in the kth state (here we take Z(k) as an n-dimensional column vector ), the calculation equation of the observation matrix is Z(k) = H*X(k) + V(k), expressed in matrix form as,
[0103] That is
[0104]
[0105] (E) The state prediction generation part calculates the transmission matrix A and the optimal estimation value Xkf(k-1) of the (k-1)th state calculated in (H) (here we take Xkf(k-1) as an n-dimensional column vector ) to get the predicted value X_pre(k) of the kth state (here we take X_pre(k) as an n-dimensional column vector ). The calculation equation of the state prediction is X_pre(k) = A*Xkf(k-1), expressed in matrix form as,
[0106] That is
[0107]
[0108] (F) The predicted covariance generation part is used to calculate the predicted state covariance P_pre(k) (here we take P_pre(k) as an n*n dimensional matrix ), the predicted state covariance is calculated by the covariance P(k-1) of the (k-1)th state calculated in (I) (here we take P(k-1) as an n*n dimensional matrix ) and the transmission matrix A calculated in (A) ) and the process noise covariance matrix Q calculated in (L) (here we take Q as an n*n dimensional matrix ), and its calculation equation is P_pre(k) = A*P(k-1)*A'+Q, expressed in matrix form as, Here then
[0109] (G) The Kalman Gain calculation section is used to calculate the Kalman Gain Kg (Kalman Gain) (here we take Kg as an n-dimensional column vector ) and the kth state prediction covariance calculated by (F) and the observation matrix calculated by (B) and the measurement noise covariance R calculated by (M) to operate. The calculation equation is Kg(k) = P_pre(k)*H' *inv(H*P_pre(k)*H' + R). It is expressed in matrix form, Here we set
[0110] where Kgd11 = ((h11*P pre11(k) +…+h1n*P pre1n(k) )*h11+…+(h11*P pre1n(k) +…+h1n*P prenn(k) )*h1n)+R, …, Kgd1n = ((h11*P pre11(k) +…+h1n*P pre1n(k) )*hn1+…+ (h11*P pre1n(k) +…+h1n*P prenn(k) )*hnn)+R, …, Kgdn1 = ((hn1*P _pre11(k) +…+ hnn*P _pren1(k) )*h11+…+(hn1*P pre1n(k) +…+hnn*P prennk) )*h1n)+R, …, Kgdnn = ((hn1*P _pre11(k) +…+hnn*P -pren1(k) )*hn+,…+(hn1*P _pre1n(k) +…+hnn* P _prenn(k) )*hnn)+R. Then
[0111] (H) The state optimal estimation value calculation section is used to calculate the state optimal estimation value The state optimal estimation value is calculated by the state prediction value calculated in (F) and the Kalman Gain calculated in (G) and the observation measurement value calculated in (D) and the observation matrix calculated in (B) is calculated. Its calculation equation is Xkf(k) = X_pre(k) + Kg(k) * (Z(k) - H * X_pre(k)).
[0112] is expressed in matrix form as,
[0113]
[0114] The (I) state covariance update part is generated by the (G) Kalman gain and the (F) state prediction covariance generated by the (B) and the (B) observation matrix calculated by the (C) The state covariance is updated by calculation Its calculation equation is P(k) = (I - Kg(k) * H) * P_pre(k). It is expressed in matrix form as,
[0115]
[0116] where P11(k) = ((1 - (Kg11(k) * h11+... + Kg1n(k) * hn1)) * P_pre11(k) +... + (-(Kg11(k) * h1n+... + Kg1n(k) * hnn)) * P -pren1 (k)),..., P1n = ((1 - (Kg11(k) * h11+... + Kg1n(k) * hn1) * P -pre1n (k) +... + (-(Kg11(k) * h1n+... + Kg1n(k) * hnn)) * P -prenn (k))),..., Pn1(k) = ((-(Kgn1(k) * h11+... + Kgnn(k) * hn1) * P _pre11 (k) +... + (1 - (Kgn1(k) * h1n+... + Kgnn(k) * hnn)) * P -pren1 (k)),..., Pnn = ((-(Kgn1(k) * h11+... + Kgnn(k) * hn1) * P -pre1n (k) +... + (1 - (Kgn1(k) * h1n+... + Kgnn(k) * hnn)) * P -prenn (k)))
[0117] (J) (K) are respectively the noise encountered in the measurement process and the noise generated in the state transition process and The covariances of these two noises are respectively the covariances R and
[0118] The actual situation and experience formula calculation are generated according to the measurement, and the specific form of the formula is different due to different systems and measurement processes.
[0119] The embodiment of the application discloses a parallel Kalman filtering system, referring to Figures 1-5 As shown in the figure, including the transmission matrix generation module, observation matrix generation module, state transition calculation module, observation measurement calculation module, state prediction generation module, prediction covariance generation module, Kalman gain calculation module, state optimal estimation value calculation module, state covariance update module, noise generation module in the measurement process, noise generation module generated in the state transition process, process noise covariance generation module and measurement noise covariance generation module. The state transition calculation module, the observation measurement calculation module, the state prediction generation module, the prediction covariance generation module and the Kalman gain calculation module are all connected with the transmission matrix generation module and the observation matrix generation module. The state optimal estimation value calculation module and the state covariance update module are both connected with the observation matrix generation module. The prediction covariance generation module, the Kalman gain calculation module, the noise generation module in the measurement process and the noise generation module generated in the state transition process are all connected with the process noise covariance generation module and the measurement noise covariance generation module. The noise generation module in the measurement process and the noise generation module generated in the state transition process are respectively connected with the observation measurement calculation module and the state transition calculation module.
[0120] (1) The transmission matrix generation module (AModule) can complete the function in step (A).
[0121] (2) The observation matrix generation module (HModule) can complete the function in step (B).
[0122] (3) The State Equation Module is implemented by a plurality of independent sub-modules (SESubModule) in parallel: State Equation Sub-Module 1 (SESubModulel) calculates X(l): X(l) = system initial state. State Equation Sub-Module 2 (SESubModule2) calculates X(2): X(2) = A*X(l) + L*W(l). State Equation Sub-Module N-l (SESubModuleN-l) calculates X(N-l): X(N-l) = A*X(N-2) + L*W(N-2). State Equation Sub-Module N (SESubModuleN) calculates X(N): X(N) = A*X(N-l) + L*W(N-l). Each State Equation Sub-Module can complete the full function X(k) = A*X(k-l) + L*W(k-l) described in step (C), so N State Equation Sub-Modules can simultaneously calculate N different system states X(N). When N < K, after N sub-modules have calculated N different system states X(l), X(2)... X(N-l), X(N), a new round of calculation begins, and N different system states X(N+l), X(N+2)... X(2N-l), X(2N) are calculated, and so on until the K states of the system X(l), X(2)... X(K-l), X(K) are all calculated. When N >= K, K of the N sub-modules can directly calculate all K states of the system X(l), X(2)... X(K-l), X(K) at one time.
[0123] (4) The Observe Equation Module (OESubModule) is implemented with a parallel set of independent sub-modules (OESubModule): Observe Equation Sub-Module 1 (OESubModulel) computes Z(l): Z(l) = H*X(l) + V(l). Observe Equation Sub-Module 2 (OESubModule2) computes Z(2): Z(2) = H*X(2) + V(2). Observe Equation Sub-Module N-l (OESubModuleN-l) computes Z(N-l): Z(N-l) = H*X(N-l) + V(N-l). Observe Equation Sub-Module N (OESubModuleN) computes Z(N): Z(N) = H*X(N) + V(N). Each Observe Equation Sub-Module is capable of performing the complete function Z(k) = H*X(k) + V(k) described in step (D) so that N Observe Equation Sub-Modules can simultaneously compute N different values of the observed equation Z(N). When N < K, after N sub-modules have computed N different values of the observed equation Z(l), Z(2)... Z(N-l), Z(N), a new round of computation begins and N different values of the observed equation Z(N+l), Z(N+2)... Z(2N-l), Z(2N) are computed, and so on until the K values of the observed equation X(l), X(2)... X(K-l), X(K) have all been computed. When N >= K, K of the N sub-modules can be used to directly compute all K values of the observed equation Z(l), Z(2)... Z(K-l), Z(K) at one time.
[0124] (5) The State Prediction Generation Module (SPESubModule) is implemented by a parallel multiple independent sub-modules (SPESubModule): State Prediction Generation Sub-Module 1 (SPESubModule1) calculates X_pre(1): X_pre(1) = System State Prediction Initial Value. State Prediction Generation Sub-Module 2 (SPESubModule2) calculates X_pre(2): X_pre(2) = A*Xkf(1). State Prediction Generation Sub-Module N-1 (SPESubModuleN-1) calculates X_pre(N-1): X_pre(N-1) = A*Xkf(N-2). State Prediction Generation Sub-Module N (SPESubModuleN) calculates X_pre(N): X_pre(N) = A*Xkf(N-1). Each State Prediction Generation Sub-Module can complete the full function X_pre(k) = A*Xkf(k-1) described in step (E), so N State Prediction Generation Sub-Modules can simultaneously calculate N different state prediction values X_pre(N). When N < K, after N sub-modules calculate N different state prediction values X_pre(1), X_pre(2)…X_pre(N-1), X_pre(N), a new round of calculation begins, and N different state prediction values X_pre(N+1), X_pre(N+2)…X_pre(2N-1), X_pre(2N) are calculated, and so on until the K state prediction values X_pre(1), X_pre(2)…X_pre(K-1), X_pre(K) of the system are calculated. When N >= K, K sub-modules of N sub-modules can directly calculate the K state prediction values X_pre(1), X_pre(2)…X_pre(K-1), X_pre(K) of the system at one time.
[0125] (6) The Convariance Predicte Equation Module is implemented with parallel multiple independent sub-modules (CPE SubModule). The SPESubModule 1 calculates P_pre(1): P_pre(1) = System Predicted Covariance Initial Value. The SPESubModule 2 calculates P_pre(2): P_pre(2) = A*P(1)*A' + Q. The SPESubModule N-1 calculates P_pre(N-1): P_pre(N-1) = A*P(N-2)*A' + Q. The SPESubModule N calculates P_pre(N): P_pre(N) = A*P(N-1)*A' + Q. Each of the SPESubModules can complete the full function P_pre(k) = A*P(k-1)*A' + Q described in step (F). So N SPESubModules can calculate N different predicted covariance values P_pre(N) simultaneously. When N < K, after N SPESubModules calculate N different predicted covariance values P_pre(1), P_pre(2)…P_pre(N-1), P_pre(N), a new round of calculation begins and N different predicted covariance values P_pre(N+1), P_pre(N+2)…P_pre(2N-1), P_pre(2N) are calculated. This continues until K predicted covariance values P_pre(1), P_pre(2)…P_pre(K-1), P_pre(K) are calculated. When N >= K, K of the N SPESubModules can calculate K predicted covariance values P_pre(1), P_pre(2)…P_pre(K-1), P_pre(K) simultaneously.
[0126] (7) The Kalman gain calculation module (Kalman Gain Equation Module) is implemented by a plurality of independent sub-modules (KGESubModule) in parallel: Kalman gain calculation sub-module 1 (KGESubModule 1) calculates Kg(1): Kg(1) = initial value of Kalman gain. Kalman gain calculation sub-module 2 (KGESubModule 2) calculates Kg(2): Kg(2) = P_pre(2)*H'*inv(H*P_pre(2)*H'+R). Kalman gain calculation sub-module N-1 (KGESubModule N-1) calculates Kg(N-1): Kg(N-1) = P_pre(N-1)*H'*inv(H*P_pre(N-1)*H'+R). Kalman gain calculation sub-module N (KGESubModule N) calculates Kg(N): Kg(N) = P_pre(N)*H'*inv(H*P_pre(N)*H'+R). Each Kalman gain calculation sub-module can complete the full function Kg(k) = P_pre(k)*H'*inv(H*P_pre(k)*H'+R) described in step (G), so that N Kalman gain calculation sub-modules can simultaneously calculate N different Kalman gain values Kg(N). When N < K, after N sub-modules calculate N different Kalman gain values Kg(1), Kg(2)…Kg(N-1), Kg(N), a new round of calculation begins, and N different Kalman gain values Kg(N+1), Kg(N+2)…Kg(2N-1), Kg(2N) are calculated, and so on until the K Kalman gain values Kg(1), Kg(2)…Kg(K-1), Kg(K) of the system are calculated. When N >= K, K sub-modules of the N sub-modules can directly calculate the K Kalman gain values Kg(1), Kg(2)…Kg(K-1), Kg(K) of the system at one time.
[0127] (8) The State Update Equation Module is implemented by a plurality of independent sub-modules (SUE Sub-Module) in parallel: State Update Equation Sub-Module 1 (KGE Sub-Module 1) calculates Xkf(1): Xkf(1) = State Update Equation Initial Value. State Update Equation Sub-Module 2 (KGE Sub-Module 2) calculates Xkf(2): Xkf(2) = X_pre(2) + Kg(2)*(Z(2) - H*X_pre(2). State Update Equation Sub-Module N-1 (KGE Sub-Module N-1) calculates Xkf(N-1): Xkf(N-1) = X_pre(N-1) + Kg(N-1)*(Z(N-1) - H*X_pre(N-1). State Update Equation Sub-Module N (KGE Sub-Module N) calculates Xkf(N): Xkf(N) = X_pre(N) + Kg(N)*(Z(N) - H*X_pre(N). Each State Update Equation Sub-Module can complete the full function Xkf(k) = X_pre(k) + Kg(k)*(Z(k) - H*X_pre(k) described in step (H), so N State Update Equation Sub-Modules can simultaneously calculate N different State Update Equations Xkf(N). When N < K, after N different State Update Equations Xkf(1), Xkf(2)... Xkf(N-1), Xkf(N) are calculated by N Sub-Modules, a new round of calculation begins, and N different State Update Equations Xkf(N+1), Xkf(N+2)... Xkf(2N-1), Xkf(2N) are calculated, and so on until the K State Update Equations Xkf(1), Xkf(2)... Xkf(K-1), Xkf(K) of the system are all calculated. When N >= K, K State Update Equations Xkf(1), Xkf(2)... Xkf(K-1), Xkf(K) of the system can be directly calculated by K of the N Sub-Modules at one time.
[0128] (9) The state covariance update module (ConvarianceUpdateEquationModule) is implemented with parallel multiple independent sub-modules (CUESubModule): state covariance update sub-module 1 (CUESubModule1) calculates P(1): P(1) = state covariance initial value. State covariance update sub-module 2 (CUESubModule2) calculates P(2): P(2) = (I - Kg(2)*H)*P_pre(2). State covariance update sub-module N-1 (CUESubModuleN-1) calculates P(N-1): P(N-1) = (I - Kg(N-1)*H)*P_pre(N-1). State covariance update sub-module N (CUESubModuleN) calculates P(N): P(N) = (I - Kg(N)*H)*P_pre(N). Each state covariance update sub-module can complete the full function P(k) = (I - Kg(k)*H)*P_pre(k) described in step (I), so N state covariance update sub-modules can simultaneously calculate N different state covariance update values P(N). When N < K, after N sub-modules calculate N different state covariance update values P(1), P(2)…P(N-1), P(N), a new round of calculation begins, and N different state covariance update values P(N+1), P(N+2)…P(2N-1), P(2N) are calculated, and so on until the K state covariance update values P(1), P(2)…P(K-1), P(K) of the system are all calculated. When N >= K, K sub-modules of N sub-modules can directly calculate the K state covariance update values P(1), P(2)…P(K-1), P(K) of the system at one time.
[0129] (10) The VGenerateModule is implemented with a plurality of independent submodules (VGSubModule) in parallel: a VGenerateSubmodule 1 (VGSubModule 1) that calculates V(1); a VGenerateSubmodule 2 (VGSubModule 2) that calculates V(2); a VGenerateSubmodule N-1 (VGSubModule N-1) that calculates V(N-1); and a VGenerateSubmodule N (VGSubModule N) that calculates V(N). Each VGenerateSubmodule is capable of performing the complete function described in step (J), so that N VGenerateSubmodules can simultaneously calculate N different V(N) values. When N < K, after N VGenerateSubmodules have calculated N different V(N) values V(1), V(2)... V(N-1), V(N), a new round of calculation begins, and N VGenerateSubmodules calculate N different V(N) values V(N+1), V(N+2)... V(2N-1), V(2N), and so on, until K V(N) values V(1), V(2)... V(K-1), V(K) have been calculated. When N >= K, K VGenerateSubmodules can directly calculate K V(N) values V(1), V(2)... V(K-1), V(K) at one time.
[0130] (11) The WGenerateModule is implemented with a plurality of independent submodules (WGSubModule) in parallel: a WSubModulel that calculates W(l); a WSubModule2 that calculates W(2); a WSubModuleN-1 that calculates W(N-l); and a WSubModuleN that calculates W(N). Each of the WSubModules can perform the complete function described in step (K), so that N WSubModules can simultaneously calculate N different W(N). When N < K, after N WSubModules have calculated N different W(l), W(2),..., W(N-l), W(N), a new round of calculation begins, and N WSubModules calculate N different W(N+l), W(N+2),..., W(2N-l), W(2N). This continues until the K W(l), W(2),..., W(K-l), W(K) of the system have all been calculated. When N >= K, K WSubModules can directly calculate the K W(l), W(2),..., W(K-l), W(K) of the system at one time.
[0131] (12) The QGenerateModule can perform the function in step (L).
[0132] (13) The RGenerateModule can perform the function in step (M).
[0133] Specifically, referring to FIG. 1, first, the module names in the figure are explained as follows: Figure 5
[0134] AM (AModule): transmission matrix generation module. HM (HModule): observation matrix generation module. SESM (State Equation SubModule): state transition calculation submodule. OESM (Observe Equation SubModule): observation measurement calculation submodule. SPESM (State Predict Equation SubModule): state prediction generation submodule. CPESM (Convariance Predict Equation SubModule): predicted covariance generation submodule. KGESM (Kalman Gain Equation SubModule): Kalman gain calculation submodule. SUESM (State Update Equation SubModule): state optimal estimation value calculation submodule. CUESM (Convariance Update Equation SubModule): state covariance update submodule. VGSM (VGenerate SubModule): measurement noise generation submodule. WGSM (WGenerate SubModule): transmission noise generation submodule. RGM (RGenerate Module): process noise covariance generation module. QGM (QGenerate Module): measurement noise covariance generation module.
[0135] The above description of each module or submodule is in terms of the function it performs. If it is an ASIC, the above each module or submodule corresponds to the corresponding functional module of the special integrated circuit; if it is a multi-core CPU / GPU, the above each module or submodule corresponds to one or more CPU / GPU cores used to complete the function; if it is an FPGA, the above each module or submodule corresponds to one or more FPGA units (such as DSP block, LUT, BRAM, etc.) that implement the function.
[0136] The data transfer and operation method among parallel Kalman filtering submodules includes the following steps:
[0137] 1) AM (AModule-Transmission Matrix Generation Module) calculates the generated transmission matrix Through paths 1_1, 1_2, …, 1_n-1, 1_n, it is transmitted to SESM (State Equation SubModule- State Transition Calculation Submodule) 1, SESM2, …, SESMn-1, SESMn.
[0138] 2) AM calculates the generated transmission matrix Through paths 2_1, 2_2,..., 2_n-1, 2_n to SPESM (State Predict Equation SubModule) 1, SPESM2,..., SPESMn-1, SPESMn respectively.
[0139] 3) AM calculates generated transmission matrix Through paths 3_1, 3_2,..., 3_n-1, 3_n to CPESM (Convariance Predict Equation SubModule) 1, CPESM2,..., CPESMn-1, CPESMn respectively.
[0140] 4) HM (HModule) calculates generated observation matrix Through paths 4_1, 4_2,..., 4_n-1, 4_n to OESM (Observe Equation SubModule) 1, OESM2,..., OESMn-1, OESMn respectively.
[0141] 5) HM calculates generated observation matrix Through paths 5_1, 5_2,..., 5_n-1, 5_n to KGESM (Kalman Gain Equation SubModule) 1, KGESM2,..., KGESMn-1, KGESMn respectively.
[0142] 6) SESM1 calculates generated X(1) = system initial state Through path 6_1 to SESM2, SESM2 calculates generated Through path 6_2 to SESM3, SESMn-1 calculates generated Through path 6_n-1 to SESMn, SESMn calculates generated Through path 6_n to SESM1 (as Figure 3 (3) When N < K, this operation is performed. When N >= K, SESMk calculates generated End operation).
[0143] 7) WGSM (QGenerate Module) 1 calculates generated Through path 7_1 to SESM2, WGSM2 calculates generated WGSMn-1 computes the generation through path 7_2 to SESM3 WGSMn computes the generation through path 7_n-1 to SESMn-1 WGSMk computes the generation through path 7_n to SESMl (When N < K, this is done. When N >= K, WGSMk computes the generation and ends the operation).
[0144] 8) SESMl computes the generation SESM2 computes the generation through path 8_1 to OESMl SESMn-1 computes the generation through path 8_2 to OESM2 SESMn computes the generation through path 8_n-1 to OESMn-1 SESMn computes the generation through path 8_n to OESMn.
[0145] 9) VGSM (VGenerateSubModule) 1 computes the generation VGSM2 computes the generation through path 9_1 to OESMl VGSMn-1 computes the generation through path 9_2 to OESM2 VGSMn computes the generation through path 9_n-1 to OESMn-1 VGSMn computes the generation through path 9_n to OESMn.
[0146] 10) SUESM (StateUpdateEquationSubModule) 1 computes the generation SUESM2 computes the generation through path 10_1 to SPESM (StatePredicteEquationSubModule) 2 SUESMn-1 computes the generation through path 10_2 to SPESM3 SUESMn computes the generation through path 10_n-1 to SPESMn SUESMk computes the generation through path 10_n to SPESMl (When N < K, this is done. When N >= K, SUESMk computes the generation and ends the operation).
[0147] 11) CUESM (ConvarianceUpdateEquationSubModule) 1 computes the generation transmitted to CPESM2 through path 11_1, CUESM2 calculates the generated transmitted to CPESM3 through path 11_2, CUESMn-1 calculates the generated transmitted to CPESMn through path 11_n-1, CUESMn calculates the generated transmitted to CPESM1 through path 11_n (when N < K, this operation is performed. When N >= K, CUESMk calculates the generated and ends the operation).
[0148] 12) QGM (QGenerate Module - Measurement Noise Covariance Generation Module) calculates the generated process noise covariance matrix transmitted to CPESM1, CPESM2,..., CPESMn-1, CPESMn through paths 12_1, 12_2,..., 12_n-1, 12_n, respectively.
[0149] 13) CPESM1 calculates the generated P_pre(1) = System Prediction Covariance Initial Value transmitted to KGESM1 through path 13_1, CPESM2 calculates the generated transmitted to KGESM2 through path 13_2, CPESMn-1 calculates the generated transmitted to KGESMn-1 through path 13_n-1, CPESMn calculates the generated transmitted to KGESMn through path 13_n.
[0150] 14) RGM (RGenerate Module - Process Noise Covariance Generation Module) calculates the generated measurement covariance R transmitted to KGESM1, KGESM2,..., KGESMn-1, KGESMn through paths 14_1, 14_2,..., 14_n-1, 14_n, respectively.
[0151] 15) SPESM1 calculates the generated transmitted to SUESM1 through path 15_1, SPESM2 calculates the generated transmitted to SUESM2 through path 15_2, SPESMn-1 calculates the generated transmitted to SUESMn-1 through path 15_n-1, SPESMn calculates the generated transmitted to SUESMn through path 15_n.
[0152] 16) KGESM1 calculates the generated Kg Kgesm2 computes the generated Kg and sends it to SUESM2 through path 16_2 Kgesmn-1 computes the generated Kg and sends it to SUESMn-1 through path 16_n-1 Kgesmn computes the generated Kg and sends it to SUESMn through path 16_n
[0153] 17) Oesml computes the generated Oesm2 computes the generated and sends it to SUESM2 through path 17_2 Oesmn-1 computes the generated and sends it to SUESMn-1 through path 17_n-1 Oesmn computes the generated and sends it to SUESMn through path 17_n
[0154] 18) Kgesml computes the generated Kg Kgesm2 computes the generated Kg and sends it to CUESM2 through path 18_2 Kgesmn-1 computes the generated Kg and sends it to CUESMn-1 through path 18_n-1 Kgesmn computes the generated Kg and sends it to CUESMn through path 18_n
[0155] 19) Cpesml computes the generated Cpesm2 computes the generated and sends it to CUESM2 through path 19_2 Cpesmn-1 computes the generated and sends it to CUESMn-1 through path 19_n-1 Cpesmn computes the generated and sends it to CUESMn through path 19_n
[0156] 20) RGM computes the generated measurement noise covariance R and sends it to VGSMl, VGSM2,..., VGSMn-1, VGSMn through paths 20_1, 20_2,..., 20_n-1, 20_n respectively.
[0157] 21) QGM computes the generated process noise covariance matrix The path 21_1, 21_2,..., 21_n-1, 21_n is transmitted to the WGS M1, WGS M2,..., WGS Mn-1, WGS Mn respectively.
[0158] In a preferred embodiment of the present application, 20000 state points are tracked by using parallel Kalman filter. Doppler frequency offset DetaTheta, frequency offset change rate Omega0, and frequency offset change rate change rate Omega1 are selected as the state quantities of the system, so the system in this example is a third-order system, and detav=1e-3, detaw=1e-3, and Ts (sampling time)=1e-3 are taken. The observation matrix generation module (HM) calculates the empirical formula of the observation matrix as H=[1 Ts / 2 Ts^2 / 6]=[1 5.0e-04 1.67e-07]. The transmission matrix generation module (AM) calculates the empirical formula of the observation matrix as The measurement noise covariance generation module (QGM) calculates the empirical formula of the process noise covariance matrix as
[0159] The process noise covariance generation module (RGM) calculates the empirical formula of the measurement noise covariance as R=detaw^2 / Ts^4=1.0e+06. The measurement noise generation submodule (VGSM) calculates the measurement noise V=N(0, R). The transmission noise generation submodule (WGSM) calculates the transmission noise W=N(0, Q). The initial value of the system state quantity is set as X(1)=[0 -20*e3 0], the initial value of the system observation quantity is set as Z(1)=[0 0 0], and the initial value of the system state covariance is set as P(1)=[1 0 0; 0 1 0; 0 0 1]. The initial value of the prediction covariance is set as P_pre(1)=[1 0 0; 0 1 0; 0 0 1]. The initial value of the state optimal estimate Xkf(1)=[0 -20*10^(3) 0], and the initial value of the state prediction X_pre(1)=[0 0 0]. Referring to Figure 6 As shown in FIG. 5, n=500 is taken, that is, each step of the Kalman filter is operated in parallel by 500 submodules.
[0160] As described above, in this example, the SESMk completes X(k)=A*X(k-1)+L*W(k-1), which is expressed in matrix form as,
[0161] which is expressed in equation form as,
[0162] x1(k) = a11 * x1(k - 1) + a12 * x2(k - 1) + a13 * x3(k - 1) + 111 * w1(k - 1) + 112 * w2(k - 1) + 113 * w3(k - 1),
[0163] x2(k) = a21 * x1(k - 1) + a22 * x2(k - 1) + a23 * x3(k - 1) + 121 * w1(k - 1) + 122 * w2(k - 1) + 123 * w3(k - 1),
[0164]
[0165] x1(k) = x1(k - 1) + 1.0e-03 * x2(k - 1) + 5.0e-07 * x3(k - 1) + w1(k - 1),
[0166] x2(k) = + x2(k - 1) + 1.0e-03 * x3(k - 1) + w2(k - 1),
[0167]
[0168] (A) As described above in this example, OESMk accomplishes Z(k) = H * X(k) + V(k), expressed in matrix form as
[0169]
[0170] expressed in equation form as,
[0171]
[0172] (B) As described above in this example, SPESMk accomplishes X_pre(k) = A * Xkf(k - 1), expressed in matrix form as expressed in equation form as,
[0173] X_pre1(k) = a11 * Xkf1(k - 1) + a12 * Xkf2(k - 1) + a13 * Xkf3(k - 1),
[0174] X_pre2(k) = a21 * Xkf1(k - 1) + a22 * Xkf2(k - 1) + a23 * Xkf3(k - 1),
[0175]
[0176] X_pre1(k) = Xkf1(k - 1) + 1.0e-03 * Xkf2(k - 1) + 5.0e-07 * Xkf3(k - 1),
[0177] X_pre2(k)=Xkf2(k-1)+1.0e-03*Xkf3(k-1),
[0178]
[0179] (C) As described above in this example, CPESMk completes the calculation equation P_pre(k)=A*P(k-1)*A'+Q, which can be expressed in matrix form as follows:
[0180]
[0181] Expressed in equation form,
[0182] P_pre11(k)=
[0183] (A11*P11(k-1)+A12*P21(k-1)+A13*P31(k-1))*A11+(A11*P12(k-1)+A12*P22(k -1)+A13*P32(k-1))*A12+(A11*P13(k-1)+A12*P23(k-1)+A13*P33(k-1))*A13+Q 11;
[0184] P_pre12(k)=
[0185] (A11*P11(k-1)+A12*P21(k-1)+A13*P31(k-1))*A21+(A11*P12(k-1)+A12*P22(k -1)+A13*P32(k-1))*A22+(A11*P13(k-1)+A12*P23(k-1)+A13*P33(k-1))*A23+Q 12;
[0186] P_pre13(k)=
[0187] (A11*P11(k-1)+A12*P21(k-1)+A13*P31(k-1))*A31+(A11*P12(k-1)+A12*P22(k -1)+A13*P32(k-1))*A32+(A11*P13(k-1)+A12*P23(k-1)+A13*P33(k-1))*A33+Q 13;
[0188] P_pre21(k)=
[0189] (A21*P11(k-1) + A22*P21(k-1) + A23*P31(k-1))*A31 + (A21*P12(k-1) + A22*P22(k-1) + A23*P32(k-1))*A32 + (A21*P13(k-1) + A22*P23(k-1) + A23*P33(k-1))*A33 + Q 23;
[0190] P_pre21(k) = (A11*P11(k-1) + A12*P21(k-1) + A13*P31(k-1))*A21 + (A11*P12(k-1) + A12*P22(k-1) + A13*P32(k-1))*A22 + (A11*P13(k-1) + A12*P23(k-1) + A13*P33(k-1))*A23 + Q 21;
[0191] (A21*P11(k-1) + A22*P21(k-1) + A23*P31(k-1))*A21 + (A21*P12(k-1) + A22*P22(k-1) + A23*P32(k-1))*A22 + (A21*P13(k-1) + A22*P23(k-1) + A23*P33(k-1))*A23 + Q 22;
[0192] P_pre23(k) = (A11*P11(k-1) + A12*P21(k-1) + A13*P31(k-1))*A31 + (A11*P12(k-1) + A12*P22(k-1) + A13*P32(k-1))*A32 + (A11*P13(k-1) + A12*P23(k-1) + A13*P33(k-1))*A33 + Q 23;
[0193] (A21*P11(k-1) + A22*P21(k-1) + A23*P31(k-1))*A31 + (A21*P12(k-1) + A22*P22(k-1) + A23*P32(k-1))*A32 + (A21*P13(k-1) + A22*P23(k-1) + A23*P33(k-1))*A33 + Q 23;
[0194] P_pre31(k) = (A11*P11(k-1) + A12*P21(k-1) + A13*P31(k-1))*A21 + (A11*P12(k-1) + A12*P22(k-1) + A13*P32(k-1))*A22 + (A11*P13(k-1) + A12*P23(k-1) + A13*P33(k-1))*A23 + Q 21;
[0195] (A31*P11(k-1) + A32*P21(k-1) + A33*P31(k-1))*A11 + (A31*P12(k-1) + A32*P22(k-1) + A33*P32(k-1))*A12 + (A31*P13(k-1) + A32*P23(k-1) + A33*P33(k-1))*A13 + Q 31;
[0196] P_pre32(k) = (A11*P11(k-1) + A12*P21(k-1) + A13*P31(k-1))*A31 + (A11*P12(k-1) + A12*P22(k-1) + A13*P32(k-1))*A32 + (A11*P13(k-1) + A12*P23(k-1) + A13*P33(k-1))*A33 + Q 32;
[0197] (A31*P11(k-1) + A32*P21(k-1) + A33*P31(k-1))*A21 + (A31*P12(k-1) + A32*P22(k-1) + A33*P32(k-1))*A22 + (A31*P13(k-1) + A32*P23(k-1) + A33*P33(k-1))*A23 + Q 32;
[0198] P_pre33(k) = (A11*P11(k-1) + A12*P21(k-1) + A13*P31(k-1))*A31 + (A11*P12(k-1) + A12*P22(k-1) + A13*P32(k-1))*A32 + (A11*P13(k-1) + A12*P23(k-1) + A13*P33(k-1))*A33 + Q 23;
[0199] (A31*P11(k-1)+A32*P21(k-1)+A33*P31(k-1))*A31+(A31*P12(k-1)+A32*P22(k-1)+A33*P32(k-1))*A32+(A31*P13(k-1)+A32*P23(k-1)+A33*P33(k-1))*A33+Q33;
[0200] (D) As described above in this example, the KGESMk completes the calculation equation as Kg(k) = P_pre(k)*H' *inv(H*P_pre(k)*H' + R). Represented in matrix form as,
[0201]
[0202] Represented in equation form, let
[0203] kgd = (P_pre11(k)*H1*H1 + P_pre21(k)*H1*H2 + P_pre31(k)*H1*H3 + P_pre12(k)*H1*H2 + P_pre22(k)*H2*H2 + P_pre32(k)*H2*H3 + P_pre13(k)*H1*H3 + P_pre23(k)*H2*H3 + P_pre33(k)*H3*H3) + R;
[0204] kgn1(k) = P_pre11(k)*H1 + P_pre12(k)*H2 + P_pre13(k)*H3;
[0205] kgn2(k) = P_pre21(k)*H1 + P_pre22(k)*H2 + P_pre23(k)*H3;
[0206] kgn3(k) = P_pre31(k)*H1 + P_pre32(k)*H2 + P_pre33(k)*H3;
[0207] Kg1 = kgn1 / kgd;
[0208] Kg2 = kgn2 / kgd;
[0209] Kg3 = kgn3 / kgd;
[0210] kgd = (P_pre11(k) + P_pre21(k)*(Ts / 2) + P_pre31(k)*(Ts^2 / 6) + P_pre12(k)*(Ts / 2) + P_pre22(k)*(Ts / 2)*(Ts / 2) + P_pre32(k)*(Ts / 2)*(Ts^2 / 6) + P_pre13(k)*(Ts^2 / 6) + P_pre23(k)*(Ts / 2)*(Ts^2 / 6) + P_pre33(k)*(Ts^2 / 6)*(Ts^2 / 6)) + (detaw^2 / Ts^4);
[0211] kgn1(k) = P_pre11(k) + P_pre12(k)*(Ts / 2) + P_pre13(k)*(Ts^2 / 6);
[0212] kgn2(k) = P_pre21(k) + P_pre22(k)*(Ts / 2) + P_pre23(k)*(Ts^2 / 6);
[0213] kgn3(k) = P_pre31(k) + P_pre32(k)*(Ts / 2) + P_pre33(k)*(Ts^2 / 6);
[0214] Kg1 = kgn1 / kgd;
[0215] Kg2 = kgn2 / kgd;
[0216] Kg3 = kgn3 / kgd;
[0217] kgd = (P_pre11(k) + P_pre21(k)*(5.0e-04) + P_pre31(k)*(1.67e-07) + P_pre12(k)*( 5.0e-04) + P_pre22(k)*(5.0e-04)*(5.0e-04) + P_pre32(k)*(5.0e-04)*(1.67e-07) + P_pre13(k)*(1.67e-07) + P_pre23(k)*(5.0e-04)*(1.67e-07) + P_pre33(k)*(1.67e-07)*(1.67e-07)) + ( 1.0e+06);
[0218] kgn1(k) = P_pre11(k) + P_pre12(k)*(5.0e-04) + P_pre13(k)*(1.67e-07);
[0219] kgn2(k) = P_pre21(k) + P_pre22(k)*(5.0e-04) + P_pre23(k)*(1.67e-07);
[0220] kgn3(k) = P_pre31(k) + P_pre32(k)*(5.0e-04) + P_pre33(k)*(1.67e-07);
[0221] Kg1 = kgn1 / kgd;
[0222] Kg2 = kgn2 / kgd;
[0223] Kg3 = kgn3 / kgd;
[0224] (E) As described above in this example, the SUESMk has its calculation equation as Xkf(k) = X_pre(k) + Kg(k)*(Z(k) - H*X_pre(k)).
[0225] In matrix form, it is expressed as,
[0226]
[0227] In equation form, it is expressed as,
[0228] Xkfl(k) = X_pre1(k) + Kg1*(Z(k) - H1*X_pre1(k) + H2*X_pre2(k) + H3*X_pre3(k));
[0229] Xkf2(k) = X_pre2(k) + Kg2*(Z(k) - H1*X_pre1(k) + H2*X_pre2(k) + H3*X_pre3(k));
[0230] Xkf3(k) = X_pre3(k) + Kg3*(Z(k) - H1*X_pre1(k) + H2*X_pre2(k) + H3*X_pre3(k));
[0231] Xkfl(k) = X_pre1(k) + Kg1*(Z(k) - X_pre1(k) + (Ts / 2)*X_pre2(k) + (Ts^2 / 6)*X_pre3(k));
[0232] Xkf2(k) = X_pre2(k) + Kg2*(Z(k) - X_pre1(k) + (Ts / 2)*X_pre2(k) + (Ts^2 / 6)*X_pre3(k));
[0233] Xkf3(k) = X_pre3(k) + Kg3*(Z(k) - X_pre1(k) + (Ts / 2)*X_pre2(k) + (Ts^2 / 6)*X_pre3(k)) ;
[0234] Xkf1(k) = X_pre1(k) + Kg1*(Z(k) - X_pre1(k) + (5.0e-04)*X_pre2(k) + (1.67e-07)*X_pre3(k)) ;
[0235] Xkf2(k) = X_pre2(k) + Kg2*(Z(k) - X_pre1(k) + (5.0e-04)*X_pre2(k) + (1.67e-07)*X_pre3(k)) ;
[0236] Xkf3(k) = X_pre3(k) + Kg3*(Z(k) - X_pre1(k) + (5.0e-04)*X_pre2(k) + (1.67e-07)*X_pre3(k)) ;
[0237] (F) As described above in this example, the CUE SMk whose calculation equation is P(k) = (I - Kg(k)*H)*P_pre(k). Represented in matrix form as,
[0238]
[0239] represented in equation form as,
[0240] P11(k) = (-1)*((Kg1*H1-1)*P_pre11(k) + Kg1*H2*P_pre21(k) + Kg1*H3*P_pre31(k)) ;
[0241] P12(k) = (-1)*((Kg1*H1-1)*P_pre12(k) + Kg1*H2*P_pre22(k) + Kg1*H3*P_pre32(k)) ;
[0242] P13(k) = (-1)*((Kg1*H1-1)*P_pre13(k) + Kg1*H2*P_pre23(k) + Kg1*H3*P_pre33(k)) ;
[0243] P21(k) = (-1)*(Kg2*H1*P_pre11(k) + (Kg2*H2-1)*P_pre21(k) + Kg2*H3*P_pre31(k)) ;
[0244] P22(k) = (-1) * (Kg2*H1*P_pre12(k) + (Kg2*H2-1) *P_pre22(k) + Kg2*H3*P_pre32(k)) ;
[0245] P23(k) = (-1) * (Kg2*H1*P_pre13(k) + (Kg2*H2-1) *P_pre23(k) + Kg2*H3*P_pre33(k)) ;
[0246] P31(k) = (-1) * (Kg3*H1*P_pre11(k) + Kg3*H2*P_pre21(k) + (Kg3*H3-1) *P_pre31(k)) ;
[0247] P32(k) = (-1) * (Kg3*H1*P_pre12(k) + Kg3*H2*P_pre22(k) + (Kg3*H3-1) *P_pre32(k)) ;
[0248] P33(k) = (-1) * (Kg3*H1*P_pre13(k) + Kg3*H2*P_pre23(k) + (Kg3*H3-1) *P_pre33(k)) ;
[0249] P11(k) = (-1) * ((Kg1-1) *P_pre11(k) + Kg1*(Ts / 2) *P_pre21(k) + Kg1*(Ts^2 / 6) *P_pre31(k)) ;
[0250] P12(k) = (-1) * ((Kg1-1) *P_pre12(k) + Kg1*(Ts / 2) *P_pre22(k) + Kg1*(Ts^2 / 6) *P_pre32(k)) ;
[0251] P13(k) = (-1) * ((Kg1-1) *P_pre13(k) + Kg1*(Ts / 2) *P_pre23(k) + Kg1*(Ts^2 / 6) *P_pre33(k)) ;
[0252] P21(k) = (-1) * (Kg2*P_pre11(k) + (Kg2*(Ts / 2)-1) *P_pre21(k) + Kg2*(Ts^2 / 6) *P_pre31(k)) ;
[0253] P22(k) = (-1) * (Kg2*P_pre12(k) + (Kg2*(Ts / 2)-1)*P_pre22(k) + Kg2*(Ts^2 / 6)*P_pre32(k)) ;
[0254] P23(k) = (-1) * (Kg2*P_pre13(k) + (Kg2*(Ts / 2)-1)*P_pre23(k) + Kg2*(Ts^2 / 6)*P_pre33(k)) ;
[0255] P31(k) = (-1) * (Kg3*P_pre11(k) + Kg3*(Ts / 2)*P_pre21(k) + (Kg3*(Ts^2 / 6)-1)*P_pre31(k)) ;
[0256] P32(k) = (-1) * (Kg3*P_pre12(k) + Kg3*(Ts / 2)*P_pre22(k) + (Kg3*(Ts^2 / 6)-1)*P_pre32(k)) ;
[0257] P33(k) = (-1) * (Kg3*P_pre13(k) + Kg3*(Ts / 2)*P_pre23(k) + (Kg3*(Ts^2 / 6)-1)*P_pre33(k)) ;
[0258] P11(k) = (-1) * ((Kg1-1)*P_pre11(k) + Kg1*(5.0e-04)*P_pre21(k) + Kg1*(1.67e-07)* P_pre31(k)) ;
[0259] P12(k) = (-1) * ((Kg1-1)*P_pre12(k) + Kg1*(5.0e-04)*P_pre22(k) + Kg1*(1.67e-07)* P_pre32(k)) ;
[0260] P13(k) = (-1) * ((Kg1-1)*P_pre13(k) + Kg1*(5.0e-04)*P_pre23(k) + Kg1*(1.67e-07)* P_pre33(k)) ;
[0261] P21(k) = (-1) * (Kg2*P_pre11(k) + (Kg2*(5.0e-04)-1)*P_pre21(k) + Kg2*(1.67e-07)* P_pre31(k)) ;
[0262] P22(k) = (-1) * (Kg2*P_pre12(k) + (Kg2*(5.0e-04) - 1) * P_pre22(k) + Kg2*(1.67e-07) * P_pre32(k)) ;
[0263] P23(k) = (-1) * (Kg2*P_pre13(k) + (Kg2*(5.0e-04) - 1) * P_pre23(k) + Kg2*(1.67e-07) * P_pre33(k)) ;
[0264] P31(k) = (-1) * (Kg3*P_pre11(k) + Kg3*(5.0e-04) * P_pre21(k) + (Kg3*(1.67e-07) - 1) * P_pre31(k)) ;
[0265] P32(k) = (-1) * (Kg3*P_pre12(k) + Kg3*(5.0e-04) * P_pre22(k) + (Kg3*(1.67e-07) - 1) * P_pre32(k)) ;
[0266] P33(k) = (-1) * (Kg3*P_pre13(k) + Kg3*(5.0e-04) * P_pre23(k) + (Kg3*(1.67e-07) - 1) * P_pre33(k)) ;
[0267] 1) Complete 6 rounds of operation, i.e. 3000 system state values are generated by the parallel Kalman filter system shown in Figure 6 .
[0268] 2) Set the initial value of the system state Omega1, i.e. X(3, 3000) = 80*10^(3), the initial value of the state optimal estimate Xkf(3, 300) = 80*10^(3), and other values remain unchanged, complete 1 round of operation, i.e. 500 system state values are generated by the parallel Kalman filter system shown in Figure 6 .
[0269] 3) Set the initial value of the system state Omega1, i.e. X(3, 3500) = 0, the initial value of the state optimal estimate Xkf(3, 3500) = 0, and other values remain unchanged, complete 11 rounds of operation, i.e. 5500 system state values are generated by the parallel Kalman filter system shown in Figure 6 .
[0270] 4) Set the initial value of the system state quantity Omega1, i.e. X(3, 9000) = -80*10^(3), the initial value of the state optimal estimation value Xkf(3, 3500) = -80*10^(3), and other values remain unchanged, complete 1 round of operation, i.e. use the parallel Kalman filter system shown in FIG. 1 to generate 500 system state values. Figure 6 .
[0271] 5) Set the initial value of the system state quantity Omega1, i.e. X(3, 9500) = 0, the initial value of the state optimal estimation value Xkf(3, 3500) = 0, and other values remain unchanged, complete 12 rounds of operation, i.e. use the parallel Kalman filter system shown in FIG. 1 to generate 6000 system state values. Figure 6 .
[0272] 6) Set the initial value of the system state quantity Omega1, i.e. X(3, 15500) = 80*10^(3), the initial value of the state optimal estimation value Xkf(3, 3500) = 80*10^(3), and other values remain unchanged, complete 1 round of operation, i.e. use the parallel Kalman filter system shown in FIG. 1 to generate 500 system state values. Figure 6 .
[0273] 7) Set the initial value of the system state quantity Omega1, i.e. X(3, 16000) = 0, the initial value of the state optimal estimation value Xkf(3, 3500) = 0, and other values remain unchanged, complete 8 rounds of operation, i.e. use the parallel Kalman filter system shown in FIG. 1 to generate 4000 system state values. Figure 6 .
[0274] At this point, we have realized the carrier tracking of 20000 states by using the parallel Kalman filter system shown in FIG. 1. Figure 6 .
[0275] Draw the Doppler frequency offset DetaTheta in the 20000 system states generated by the parallel Kalman filter system shown in FIG. 1, i.e. obtain Figure 6 . Figure 7
[0276] Draw the change rate Omega0 of the frequency offset in the 20000 system states generated by the parallel Kalman filter system shown in FIG. 1, i.e. obtain Figure 6 . Figure 8
[0277] Draw the change rate Omega1 of the change rate of the frequency offset in the 20000 system states generated by the parallel Kalman filter system shown in FIG. 1, i.e. obtain Figure 6 . Figure 9
[0278] The measurement error messure_err_x(k) = Z(k) - X(l,k) of 20000 system states and the Kalman estimation error kalman_err_x(k) = Xkf(l,k) - X(l,k) are plotted, i.e. Figure 10 .
[0279] The above-described embodiments are only preferred embodiments of the present application for fully illustrating the present application, and the protection scope of the present application is not limited thereto. Any equivalent substitutions or changes made by those skilled in the art based on the present application are within the protection scope of the present application. The protection scope of the present application is subject to the claims.
Claims
1. A parallel Kalman filtering method, characterized by: The Kalman operation process is divided into multiple steps, and each step is operated in parallel while each step is operated simultaneously by multiple parallel sub-modules; The Kalman operation process is at least divided into the following steps: S1, generating a transmission matrix A; S2, generating an observation matrix H; S3, calculating the state transition of the system; S4, calculating the observation measurement; S5, generating a state prediction; S6, generating a prediction covariance P_pre(k); S7, calculating a Kalman gain; S8, calculating a state optimal estimate; S9, updating the state covariance; S10, generating noise V encountered in the measurement process and noise W generated in the state transition process; The Kalman operation process is at least divided into the following steps: S1, generating a transmission matrix A; S2, generating an observation matrix H; S3, calculating the state transition of the system; S4, calculating the observation measurement; S5, generating a state prediction; S6, generating a prediction covariance P_pre(k); S7, calculating a Kalman gain; S8, calculating a state optimal estimate; S9, updating the state covariance; S10, generating noise V encountered in the measurement process and noise W generated in the state transition process; (A) Calculate the transmission matrix A: Take A as an n*n dimensional matrix (B) Computing the observation matrix H: Take H as an n*n matrix (C) calculating the state value of the system: the state transition of the system is calculated by the state transition equation X(k) = A * X(k-1) + L * W(k-1); L is a state transition noise adjustment parameter matrix; (D) Computing the observed measurement: using the matrix calculated in step (B) and the kth state matrix multiplied and added to the measurement noise matrix V(k) in the kth state to obtain the measurement value Z(k) in the kth state, and the calculation equation of the observation matrix is: measurement value Z(k) = H * X(k) + V(k), H is the observation matrix; (E) generating a state prediction: performing matrix multiplication of the transmission matrix calculated in step (A) and the optimal estimate of the (k-1)th state calculated in step (H) to obtain a predicted value of the kth state, X_pre(k) taking Xkf(k-1) as an n-dimensional column vector taking X_pre(k) as an n-dimensional column vector the calculation equation of the state prediction is: X_pre(k) = A * Xkf(k-1) (F) generating a predicted covariance P_pre(k): taking P_pre(k) as an n*n dimensional matrix The calculation equation is P_pre(k) = A*P(k-1)*A' + Q; (G) Calculate Kalman gain: take Kalman gain Kg as n-dimensional column vector The calculation equation is Kg(k) = P_pre(k) * H' * inv(H * P_pre(k) * H' + R); H' represents the transpose matrix of the observation matrix H; (H) Calculate the state optimal estimate: take the state optimal estimate The calculation equation is Xkf(k) = X_pre(k) + Kg(k) * (Z(k) - H * X_pre(k)). (I) updating the state covariance: its calculation equation is P(k) = (I-Kg(k) * H) * P_pre(k); I represents an identity matrix; (J) (K) are the noise encountered in the measurement process and the noise generated in the state transition process, respectively and The covariances of these two noises are R and 2. A parallel Kalman filtering system, characterized by: The parallel Kalman filtering method of claim 1 comprises a transmission matrix generation module, an observation matrix generation module, a state transition calculation module, an observation measurement calculation module, a state prediction generation module, a prediction covariance generation module, a Kalman gain calculation module, a state optimal estimate calculation module, a state covariance update module, a measurement process noise generation module, a state transition process noise generation module, a process noise covariance generation module, and a measurement noise covariance generation module; the state transition calculation module, the observation measurement calculation module, the state prediction generation module, the prediction covariance generation module, and the Kalman gain calculation module are connected with the transmission matrix generation module and the observation matrix generation module; the state optimal estimate calculation module and the state covariance update module are connected with the observation matrix generation module; the prediction covariance generation module, the Kalman gain calculation module, the measurement process noise generation module, and the state transition process noise generation module are connected with the process noise covariance generation module and the measurement noise covariance generation module; the measurement process noise generation module and the state transition process noise generation module are connected with the observation measurement calculation module and the state transition calculation module, respectively.
3. The parallel Kalman filtering system of claim 2, wherein: (1) the transmission matrix generation module can complete the function in step (A); (2) the observation matrix generation module can complete the function in step (B); (3) the state transition calculation module includes multiple independent sub-modules in parallel, and each state transition calculation sub-module can complete the complete function X(k) = A * X(k-1) + L * W(k-1) in step (C); (4) the observation measurement calculation module includes multiple independent sub-modules in parallel, and each observation measurement calculation sub-module can complete the complete function Z(k) = H * X(k) + V(k) in step (D); (5) The state prediction generation module includes a plurality of independent sub-modules in parallel, each of which can complete the complete function X_pre(k) = A * Xkf(k-1) in step (E); (6) The prediction covariance generation module includes a plurality of independent sub-modules in parallel, each of which can complete the complete function P_pre(k) = A * P(k-1) * A' + Q in step (F); A' represents the transpose of matrix A; (7) The Kalman gain calculation module includes a plurality of independent sub-modules in parallel, each of which can complete the complete function Kg(k) = P_pre(k) * H' * inv(H * P_pre(k) * H' + R) in step (G); (8) The state optimal estimate calculation module includes a plurality of independent sub-modules in parallel, each of which can complete the complete function Xkf(k) = X_pre(k) + Kg(k) * (Z(k) - H * X_pre(k) in step (H); (9) The state covariance update module includes a plurality of independent sub-modules in parallel, each of which can complete the complete function P(k) = (I - Kg(k) * H) * P_pre(k) in step (I); (10) The measurement process noise generation module includes a plurality of independent sub-modules in parallel, each of which can complete the complete function in step (J); (11) The state transition process noise generation module includes a plurality of independent sub-modules in parallel, each of which can complete the complete function in step (K).
4. The parallel Kalman filtering system of claim 3, wherein: Each module or sub-module is an ASIC functional module, one or more CPU / GPU cores, or one or more FPGA units.
5. A parallel Kalman filter sub-module data transfer and operation method, characterized in that: The steps include: 1) AM calculates and generates the transmission matrix The path 1_1, 1_2, …, 1_n-1, 1_n is transmitted to the SESM1, SESM2, …, SESMn-1, SESMn respectively; AM represents the transmission matrix generation module; 2) AM computation generated transmission matrix transmitted to SPESM1, SPESM2,..., SPESMn-1, SPESMn through paths 2_1, 2_2,..., 2_n-1, 2_n respectively; 3) AM computation generated transmission matrix Through paths 3_1, 3_2, …, 3_n-1, 3_n to CPESM1, CPESM2, …, CPESMn-1, CPESMn respectively; 4) HM calculates generated observation matrix Through paths 4_1, 4_2, …, 4_n-1, 4_n, to OESM1, OESM2, …, OESMn-1, OESMn respectively; HM represents an observation matrix generation module; 5) HM computes generated observation matrix Through paths 5_1, 5_2, …, 5_n-1, 5_n to KGESM1, KGESM2, …, KGESMn-1, KGESMn respectively; 6) SESM1 calculates generated X(1) = system initial state is transmitted to SESM2 through path 6_1, SESM2 calculates generated X(2) = system state after first time step is transmitted to SESM3 through path 6_2, SESMn-1 calculates generated X(n-1) = system state after (n-1)th time step is transmitted to SESMn through path 6_n-1, SESMn calculates generated X(n) = system state after nth time step is transmitted to SESM1 through path 6_n, when N < K, this operation is performed; when N >= K, SESMk calculates generated X(K) = system state after Kth time step, and then end operation; SESM represents a state transition calculation sub-module; 7) The result generated by WGSM1 is transmitted to SESM2 via path 7_1. The result generated by WGSM2 is transmitted to SESM3 via path 7_2. The result generated by WGSMn-1 is transmitted to SESMn via path 7_n-1. The result generated by WGSMn is transmitted to SESM1 via path 7_n. This operation is performed when N < K; when N >= K, the operation ends after the result is generated by WGSMk; W(k) represents the noise vector in the k-th state transition process. 8) SESM1 computes the birth SESM2 computes the birth via path 8_1 to OESM1 SESMn-1 computes the birth via path 8_2 to OESM2 SESMn computes the birth via path 8_n-1 to OESMn-1 SESMn computes the birth via path 8_n to OESMn 9) VGSM1 computed generated VGSM2 computed generated is transmitted to OESM1 through path 9_1 VGSMn-1 computed generated is transmitted to OESM2 through path 9_2 VGSMn computed generated is transmitted to OESMn-1 through path 9_n-1 VGSMn computed generated is transmitted to OESMn through path 9_n; V(k) represents the noise vector in the kth measurement process; 10) SUESM1 calculates the generated SUESM2 calculates the generated SUESMn-1 calculates the generated SUESMn calculates the generated SUESM1, when N < K; when N >= K, SUESMk calculates the generated Xkf(k) represents the optimal estimate of the kth state; 11) CUESM1 computes generated transmitted to CPESM2 through path 11_1, CUESM2 computes generated transmitted to CPESM3 through path 11_2, CUESMn-1 computes generated transmitted to CPESMn through path 11_n-1, CUESMn computes generated transmitted to CPESM1 through path 11_n, when N < K, this operation is performed; when N >= K, CUESMk computes generated transmitted to CPESM1 through path 11_n, when N < K, this operation is performed; when N >= K, CUESMk computes generated 12) QGM computes generated process noise covariance matrix transmitted to CPESM1, CPESM2,..., CPESMn-1, CPESMn through paths 12_1, 12_2,..., 12_n-1, 12_n, respectively; 13) CPESM1 calculates generated P_pre(1) = system prediction covariance initial value and passes it to KGESM1 through path 13_1, CPESM2 calculates generated P_pre(2) and passes it to KGESM2 through path 13_2, CPESMn-1 calculates generated P_pre(n-1) and passes it to KGESMn-1 through path 13_n-1, CPESMn calculates generated P_pre(n) and passes it to KGESMn through path 13_n; 14) The measurement covariance R calculated and generated by the RGM is transmitted to the KGESM1, KGESM2, …, KGESMn-1, KGESMn through paths 14_1, 14_2, …, 14_n-1, 14_n, respectively; 15) SPESM1 computed generated transmitted to SUESM1 through path 15_1, SPESM2 computed generated transmitted to SUESM2 through path 15_2, SPESMn-1 computed generated transmitted to SUESMn-1 through path 15_n-1, SPESMn computed generated transmitted to SUESMn through path 15_n; 16) KGESM1 computed generated transmitted to SUESM1 via path 16_1, KGESM2 computed generated transmitted to SUESM2 via path 16_2, KGESMn-1 computed generated transmitted to SUESMn-1 via path 16_n-1, KGESMn computed generated transmitted to SUESMn via path 16_n; 17) OESM1 computed generated OESM2 computed generated is communicated to SUESM1 through path 17_1 OESMn-1 computed generated is communicated to SUESM2 through path 17_2 OESMn computed generated is communicated to SUESMn-1 through path 17_n-1 OESMn-1 computed generated is communicated to SUESMn through path 17_n 18) KGESM1 computed generated transmitted to CUESM1 via path 18_1, KGESM2 computed generated transmitted to CUESM2 via path 18_2, KGESMn-1 computed generated transmitted to CUESMn-1 via path 18_n-1, KGESMn computed generated transmitted to CUESMn via path 18_n; 19) CPESM1 computes generated transmitted to CUESM1 through path 19_1, CPESM2 computes generated transmitted to CUESM2 through path 19_2, CPESMn-1 computes generated transmitted to CUESMn-1 through path 19_n-1, CPESMn computes generated transmitted to CUESMn through path 19_n; 20) The measurement noise covariance R calculated and generated by the RGM is transmitted to the VGSM1, VGSM2, …, VGSMn-1, VGSMn through paths 20_1, 20_2, …, 20_n-1, 20_n, respectively; 21) QGM calculates generated process noise covariance matrix transmitted to the WGS M1, WGS M2,..., WGS Mn-1, WGS Mn through paths 21_1, 21_2,..., 21_n-1, 21_n, respectively.
6. The parallel Kalman filter sub-module data transfer and operation method of claim 5, wherein: The number of state points of carrier tracking is 20000, Doppler frequency offset DetaTheta, frequency offset change rate Omega0, and frequency offset change rate change rate Omega1 are selected as the state quantities of the system, the system is a third-order system, detav=1e-3, detaw=1e-3, Ts=1e-3, The empirical formula for calculating the observation matrix of the HM is H=[1 Ts / 2 Ts^2 / 6]=[15.0e-04 1.67e-07]; the empirical formula for calculating the observation matrix of the AM is The empirical formula for calculating the process noise covariance matrix of the QGM is The empirical formula for RGM to compute the measurement noise covariance is R = detaw^2 / Ts^4 = 1.0e+06; VGSM computes the measurement noise V = N(0, R); WGSM computes the transmission noise W = N(0, Q); The initial value of the system state quantity is set as X(l) = [0 -20*e3 0], the initial value of the system observation quantity is Z(l) = [0 0 0], and the initial value of the system state covariance is The initial value of the prediction covariance is The initial value of the state optimal estimation value is Xkf(l) = [0 -20*10^(3) 0], and the initial value of the state prediction value is X_pre(l) = [0 0 0].
7. The parallel Kalman filter sub-module data transfer and operation method of claim 6, wherein: n is 500, then the SESMk completes X(k) = A * X(k-1) + L * W(k-1); (A) The OESMk completes Z(k) = H * X(k) + V(k); (C) The CPESMk completes P_pre(k) = A * P(k-1) * A' + Q; (D) The KGESMk completes Kg(k) = P_pre(k) * H' * inv(H * P_pre(k) * H' + R); (E) The SUESMk completes Xkf(k) = X_pre(k) + Kg(k) * (Z(k) - H * X_pre(k)); (F) CUESMk completes P(k) = (I - Kg(k) * H) * P_pre(k).
8. The parallel Kalman filter sub-module data passing and operating method of claim 7, wherein: comprising the steps of 1) 6 rounds of operation are completed to generate 3000 system state values; 2) the initial value of the system state quantity Omega1, i.e. X(3, 3000) = 80*10^(3), the initial value of the state optimal estimation value Xkf(3, 300) = 80*10^(3), and other values remain unchanged, are set, 1 round of operation is completed to generate 500 system state values; 3) the initial value of the system state quantity Omega1, i.e. X(3, 3500) = 0, the initial value of the state optimal estimation value Xkf(3, 3500) = 0, and other values remain unchanged, are set, 11 rounds of operation are completed to generate 5500 system state values; 4) the initial value of the system state quantity Omega1, i.e. X(3, 9000) = -80*10^(3), the initial value of the state optimal estimation value Xkf(3, 3500) = -80*10^(3), and other values remain unchanged, are set, 1 round of operation is completed to generate 500 system state values; 5) the initial value of the system state quantity Omega1, i.e. X(3, 9500) = 0, the initial value of the state optimal estimation value Xkf(3, 3500) = 0, and other values remain unchanged, are set, 12 rounds of operation are completed to generate 6000 system state values; 6) the initial value of the system state quantity Omega1, i.e. X(3, 15500) = 80*10^(3), the initial value of the state optimal estimation value Xkf(3, 3500) = 80*10^(3), and other values remain unchanged, are set, 1 round of operation is completed to generate 500 system state values; 7) the initial value of the system state quantity Omega1, i.e. X(3, 16000) = 0, the initial value of the state optimal estimation value Xkf(3, 3500) = 0, and other values remain unchanged, are set, 8 rounds of operation are completed to generate 4000 system state values.
Citation Information
Patent Citations
FPGA (Field Programmable Gata Array)-based unscented kalman filter system and parallel implementation method
CN101777887A
A self-adaptive Kalman filtering algorithm
CN109902568A
Gain parallel calculating device for kalman filter
JP1991188552A