Distributed multi-model adaptive fusion estimation method under unknown noise covariance and measurement loss probability

Abstract

The invention discloses a distributed multi-model adaptive fusion estimation method under unknown noise covariance and measurement loss probability, which comprises the following steps of: introducing a binary judgment factor obeying Bernoulli distribution into a sensor node measurement equation in consideration of a random loss behavior measured by a sensor; the PDF of the measurement likelihood of the sensor node can be expressed in the form of an index product; in the KLA-IMM framework, conjugate priori is distributed for unknown noise covariance and binary judgment factors, and a VB technology is adopted to jointly update state, noise covariance and binary judgment factor model conditions for posteriori PDF; then, in a mixing and output stage, a KLA method is adopted to fuse model conditions for estimation, and a VB-based centralized self-adaptive interactive multi-model filter is obtained; and expanding and applying the obtained centralized self-adaptive interactive multi-model filter based on the VB to a distributed scene by using a diffusion strategy based on a covariance cross technology, and then obtaining a distributed self-adaptive interactive multi-model filter based on the VB.

Description

Technical Field

[0001] The present invention belongs to the field of multi-sensor information fusion, and particularly relates to an estimation method for distributed multi-model adaptive fusion under unknown noise covariance and measurement loss probability. Background Art

[0002] The problems of state estimation and fusion widely exist in military and civilian fields such as target tracking, battlefield surveillance, space situation awareness, and intelligent transportation. For example, in swarm warfare, a UAV cluster can be regarded as a self-organizing network topology structure, and each UAV can be regarded as a network node in the network topology structure. When monitoring or attacking dynamic targets of interest on the battlefield, the information fusion technology is used to perform accurate distributed estimation of the dynamic targets of interest based on the collected information.

[0003] Filtering and smoothing techniques are the key to state estimation and fusion, and can estimate the hidden dynamic system state from the measurement information with random noise. The typical Kalman filter can obtain the optimal state estimation in a linear system based on the Gaussian noise assumption, and has been widely extended and applied to multi-sensor fusion estimation, thereby obtaining a wider sensing range and better estimation performance than a single sensor. However, the assumptions of the typical Kalman filter are too ideal. Considering the actual engineering applications, due to the complexity of the battlefield environment, the state estimation and fusion of the dynamic system may face non-ideal situations such as uncertain system models, unknown and time-varying covariance matrices, and measurement loss probability.

[0004] At present, to solve the state estimation problem of a jump Markov system under measurement random losses, Youn et al. described the measurement random loss behavior by using a binary random variable subject to a Bernoulli distribution, assigned a conjugate prior distribution to it, and then within the framework of the Interacting Multiple Model (IMM), used the Variational Bayesian (VB) technique to jointly infer the state variable, the binary random variable, and the unknown measurement loss probability, thereby deriving a VB-based adaptive interacting multiple model Kalman filter. Based on this idea, Li et al. proposed a VB-based distributed fusion estimation algorithm by means of a hybrid consensus strategy and the IMM framework based on the Kullback-Leibler average (KLA) (KLA-IMM). However, this algorithm requires an accurately known noise covariance matrix, and this assumption is too strong. Therefore, when the given covariance matrix is inaccurate, the estimation accuracy of the distributed fusion estimation algorithm proposed by Li et al. will be greatly affected. In addition, the hybrid consensus strategy requires multi-step communication and consensus iterative calculations, so it suffers from a serious computational burden. Compared with the consensus strategy, the diffusion strategy can fuse with the input measurement values in real time, with a faster convergence speed and a smaller mean square error. Summary of the Invention

[0005] Aiming at the deficiencies of the prior art, the present invention provides a multi-sensor fusion estimation method for non-ideal situations such as uncertain system models, unknown and time-varying covariance matrices, and measurement loss probabilities. And according to the method of the present invention, a novel VB-based centralized adaptive interacting multiple model filter (VB based centralized adaptive IMM filter, VB-CAIMMF) is derived. Then, by using the diffusion strategy based on the covariance intersection technique, the proposed centralized adaptive interacting multiple model filter is extended and applied to the distributed scenario, and a VB-based distributed adaptive interacting multiple model filter (VB based distributed adaptive IMM filter, VB-DAIMMF) is proposed.

[0006] To achieve the above object, the present invention adopts the following technical solutions:

[0007] An estimation method for distributed multi-model adaptive fusion under unknown noise covariance and measurement loss probability, characterized in that the estimation method includes:

[0008] Considering the random loss behavior of sensor measurements, a binary judgment factor following a Bernoulli distribution is introduced into the sensor node measurement equation, and then the PDF of the sensor node measurement likelihood can be expressed in the form of an exponential product;

[0009] Within the KLA-IMM framework, by assigning conjugate priors to the unknown noise covariance and the binary judgment factor, the VB technique is used to jointly update the state, noise covariance, and the conditional posterior PDF of the binary judgment factor model;

[0010] Then, in the mixing and output stages, the KLA method is adopted to fuse the model conditions for estimation, and a centralized adaptive interacting multiple model filter based on VB is obtained;

[0011] Furthermore, by using a diffusion strategy based on the covariance intersection technique, the obtained centralized adaptive interacting multiple model filter based on VB is extended and applied to the distributed scenario, and a distributed adaptive interacting multiple model filter based on VB is obtained.

[0012] It should be noted that the problem to be solved by the present invention can be described by a jump Markov linear system including N s sensors, that is:

[0013] x k = F k (m k )x k-1 + w k (1)

[0014]

[0015] where, represents the state vector at time k, m k represents the system model at time k, which is controlled by a homogeneous Markov chain and takes values from a finite model set M = {1, 2,..., M} according to the known transition probability matrix Λ = [p ij M×M where, p ij = P(m k = j|m k-1 = i) represents the transition probability from model i to model j, satisfying represents the state transition function corresponding to model m k at time k, represents zero-mean Gaussian noise with an unknown and time-varying process noise covariance matrix Q k , represents the measurement vector of sensor node s at time k, represents the measurement function of sensor node s corresponding to model m k at time k,​ respectively represent that the sensor node s has an unknown measurement noise covariance matrix of zero-mean Gaussian noise. Assume that the initial state x0 satisfies and x0, w k and are uncorrelated. For the sake of convenient writing, F k (m k ) and are briefly denoted as and is a binary decision factor subject to the Bernoulli distribution, then Equation (2) can be reformulated as

[0016]

[0017] It can be seen from Equation (3) that when , it means that there is no information loss in the measurement data collected by the sensor node s at time k; when , it means that the sensor node s only collects pure noise data at time k, that is, useful state information is lost.

[0018] It should be noted that obtaining the centralized adaptive interacting multiple model filter based on variational Bayesian includes the following steps:

[0019] Step 1: Based on the above Equations (1)-(3), construct a hierarchical Gaussian state space model, that is:

[0020] According to Equation (1), the state transition PDF corresponding to model j can be given, that is

[0021]

[0022] In addition, for the sake of convenient derivation, Equation (3) can be reformulated into the following exponential product form:

[0023]

[0024] Since the binary decision factor obeys the Bernoulli distribution, its probability mass function can be written in the following form:

[0025]

[0026] where is an unknown and time-varying mixing coefficient, which represents the probability that there is no information loss in the measurement data received by the sensor node s at time k,

[0027] Considering that using conjugate distributions within the VB framework can greatly simplify the analysis, the beta distribution is selected to describe the unknown and time-varying mixing coefficient That is

[0028]

[0029] wherein and are shape parameters

[0030] For the unknown and time-varying process noise covariance matrix Q k , the unknown state covariance matrix is modeled as a mixture inverse Wishart distribution as follows:

[0031]

[0032] where L is the number of mixture inverse Wishart distributions and are the degrees of freedom parameter and the scale matrix corresponding to the l-th inverse Wishart distribution of model j, respectively, and χ k = [χ k,1 ,..., χ k,l ,..., χ k,L is the mixture coefficient vector, satisfying

[0033] Introduce an L-dimensional binary random variable where satisfies Then equation (8) can be expressed in the following hierarchical product form

[0034]

[0035] wherein

[0036] For the unknown and time-varying measurement noise covariance matrix , an inverse Wishart distribution is used to model it, that is:

[0037]

[0038] wherein and represent the degrees of freedom parameter and the scale matrix of sensor node s, respectively. Considering the slow time-varying characteristics of the mixture coefficient and the noise covariance matrix, a heuristic method can be used to describe and , that is:

[0039]

[0040] where ρ ∈ (0, 1] is the forgetting factor is the adjustment parameter The l-th nominal state noise covariance matrix, and is the unknown process noise covariance matrix Q corresponding to model j at time k - 1 k , the mixing coefficient χ k and the unknown measurement noise covariance matrix of sensor node s and the mixing coefficient parameters related to the posterior PDF,

[0041] Let x k and x k-1 be the global latent variables, be the common model latent variable and parameter set of all sensor nodes in the sensor network, be the local latent variable and parameter set of each sensor node. The IMM filter obtains the conditional posterior PDF p(Ξ k |m k = j, Z 1:k ) corresponding to different models by running a set of sub - filters in parallel, and obtains the final posterior PDF p(Ξ k |m k = j, Z 1:k ) by further fusing the conditional posterior PDFs p(Ξ k |m 1:k ) of different models, where According to Equations (4) - (11), the joint PDF p(Ξ k , Z 1:k |m k = j) corresponding to model j can be expressed in the following form:

[0042]

[0043] where, and are the initial state mean vector and the corresponding error covariance matrix corresponding to model j at time k - 1,

[0044] Step 2: Model - conditional initialization

[0045] Assume that the posterior PDF corresponding to model i at time k - 1 has been obtained, that is

[0046]

[0047] According to Bayes' rule and the law of total probability, the mixing probability can be calculated as follows:

[0048]

[0049] Similarly, the mixture PDF p(x k-1 |m k =j,Z 1:k-1 ) can be calculated as follows:

[0050]

[0051] Similarly, the mixture PDF p(χ k-1 |m k =j,Z 1:k-1 ), p(Q k-1 |m k =j,Z 1:k-1 ) and are as follows:

[0052]

[0053] Using the KLA fusion mode, the sum terms on the right - hand side of Eqs. (15) - (19) can be approximated as a single term, i.e.,

[0054]

[0055]

[0056] where and can be calculated in the following way:

[0057]

[0058] Step 3: Model conditional filtering

[0059] Since the variables in Ξ k are coupled with each other, there is no analytical solution for the conditional posterior PDF p(Ξ k |m k =j,Z 1:k ). To solve this problem, the VB technique is adopted. By choosing a distribution Q j (Ξ k ) corresponding to model j to approximate the conditional posterior PDF p(Ξ k |m k =j,Z 1:k ), based on the mean - field assumption, Q j (Ξ k ) can be expressed as:

[0060]

[0061] The coordinate ascent method is used to obtain the optimal solution of the variational distribution q j (Δ), i.e.,

[0062]

[0063] Among them

[0064] At the (t+1)-th (t = 0, 1,..., N-1) iteration, the common model variables and parameters The local variables and parameters of a single sensor node s And the update process of the common hidden state {x k , x k-1} is as follows:

[0065] A. Common model variables and parameters Update:

[0066] 4) Let Δ = Q k , substitute Equation (12) into Equation (27) for calculation, then there is

[0067]

[0068] Among them

[0069]

[0070] Among them

[0071]

[0072] 5) Let Substitute Equation (12) into Equation (27) for calculation, then there is

[0073]

[0074] Among them

[0075]

[0076] 6) Let Δ = χ k , substitute Equation (12) into Equation (27) for calculation, then there is

[0077]

[0078] Among them

[0079]

[0080] B. The local variables and parameters of a single sensor node s (s = 1, 2,..., N s ) Update:

[0081] 4) Let Substitute Equation (12) into Equation (27) for calculation, then there is

[0082] Among them, is a normalization constant, and can be obtained through the following formula

[0083]

[0084] 5) Let Substitute Equation (12) into Equation (27) for calculation, then there is

[0085]

[0086] Among them

[0087]

[0088] 6) Let Substitute Equation (12) into Equation (27) for calculation, then there is

[0089]

[0090] Among them

[0091]

[0092] C. Update of common hidden state {x k , x k-1}:

[0093] 3) Let Δ = x k-1 Substitute Equation (12) into Equation (27) for calculation, then there is

[0094]

[0095] Among them

[0096]

[0097] Among them

[0098]

[0099] 4) Let Δ = x k Substitute Equation (12) into Equation (27) for calculation, then there is

[0100]

[0101] Among them

[0102]

[0103] Among them, and is calculated as follows

[0104]

[0105] where

[0106]

[0107] In the (t + 1)-th (t = 0, 1,..., N - 1) iteration, the required mathematical expectation can be obtained in the following way, and can be calculated as follows:

[0108]

[0109] Considering q j,(t+1) (Q k ) and is an inverse Wishart distribution, it can be known that and are Wishart distributions, that is

[0110]

[0111]

[0112] Then E j,(t+1) [log|Q k |] and can be calculated as

[0113]

[0114] In addition, according to equations (31) - (38), the required mathematical expectations E j,(t+1) [χ k , E j,(t+1) [logχ k,l , and can be calculated as

[0115]

[0116] Step 4: Model probability update

[0117] Based on the Bayesian criterion, the model probability corresponding to model j can be updated as:

[0118]

[0119] where

[0120]

[0121] Among them

[0122]

[0123] Step Five: Estimation Fusion

[0124] By using the KLA method to fuse the posterior PDF of the model updated state vector, that is Then and P k∣k Can be calculated as

[0125]

[0126] Step Six: Initial Value Selection

[0127] To execute the proposed VB-CAIMMF algorithm, the initial values of the parameters need to be selected in advance. Assume that the initial mixing coefficient corresponding to model i Follows a beta distribution, that is

[0128]

[0129] It can be seen from Equation (64) that appropriate shape parameters need to be selected and In order to obtain the prior information of the initial mixing coefficient corresponding to model i For model i, the shape parameters and Can be obtained through The mathematical expectation of, that is:

[0130]

[0131] Among them, Is the probability of accurate reception of the initial nominal measurement.

[0132] Assume that the initial measurement noise covariance matrix corresponding to model i Follows an inverse Wishart distribution, that is

[0133]

[0134] According to Equation (66), for model i, Follows a Wishart distribution, that is

[0135]

[0136] For model i, the degrees of freedom parameter And the scale matrix Can be obtained through The mathematical expectation of, that is

[0137]

[0138] Among them, is the initial nominal measurement noise covariance matrix of the sensor node s.

[0139] Similar to Equations (66)-(68), the degrees of freedom parameter and the scale matrix of the initial process noise covariance matrix Q0 corresponding to model i can be obtained through the mathematical expectation of (Q0) -1 , that is

[0140]

[0141] Among them, is the first-term nominal process noise covariance matrix at the initial moment.

[0142] Assume that the initial mixing coefficient vector χ0 = [χ 0,1 ,..., χ 0,l ,..., χ 0,L corresponding to model i follows a Dirichlet distribution, that is

[0143]

[0144] Similarly, for model i, the initial prior concentration parameter can be obtained through the mathematical expectation of the initial mixing coefficient vector χ0, that is

[0145]

[0146] Among them, is expressed as the initial nominal mixing coefficient vector, and is the credibility of the l-th initial nominal process noise covariance matrix .

[0147] It should be noted that extending the centralized adaptive interacting multiple model filter to a distributed scenario includes the following steps:

[0148] Since there is no fusion center in the distributed scenario, the local latent variables and local model parameters of each sensor node will be used to replace the global latent variables and common model parameters Similar to the centralized adaptive interacting multiple model filter, the local latent variables and local model parameters can be modeled in the following form:

[0149]

[0150] Step 1: Model condition initialization

[0151] In a distributed scenario, a single sensor node s ∈ V can obtain the conditional posterior PDFs of different models in the same way as VB-CAIMMF, that is

[0152]

[0153] Step 2: Model conditional filtering

[0154] In the model conditional filtering stage of the distributed algorithm, the VB technique is used to jointly estimate the local hidden state, local variables, and parameters. For model j ∈ M, a single sensor node s In the (t + 1)-th (t = 0, 1,..., N - 1) iteration process, it mainly includes three parts, namely: local hidden state Local variable and parameter update, local hidden state Update, and CI-based diffusion fusion estimation

[0155] A. Local state Local variable and parameter update

[0156] For sensor node s, the posterior Of its local hidden state Local variables and parameters Can be updated in a similar way to equations (28) - (43), that is:

[0157] Posterior PDF of the local process noise covariance For model j ∈ M Can be updated to

[0158]

[0159] Where

[0160]

[0161] Where

[0162]

[0163] Posterior PDF of the local L-dimensional binary random variable For model j ∈ M Can be updated to

[0164]

[0165] Where

[0166]

[0167] Local L-dimensional mixing coefficient Posterior PDF: For model j ∈ M, can be updated to

[0168]

[0169] where

[0170]

[0171] local binary decision factor Posterior PDF: For model j ∈ M, can be updated to

[0172]

[0173] where, is the normalization constant, and can be obtained by the following formula

[0174]

[0175] local mixing coefficient Posterior PDF: For model j ∈ M, can be updated to

[0176]

[0177] where

[0178]

[0179] local measurement noise covariance Posterior PDF: For model j ∈ M, can be updated to

[0180]

[0181] where

[0182]

[0183] local variable Posterior PDF: For model j ∈ M, can be updated to

[0184]

[0185] where

[0186]

[0187] where

[0188]

[0189] B. Local hidden state Update

[0190] At the (t + 1)-th iteration, a single sensor node s By using its own local prior information for And the local measurement information of its neighbors for Update the posterior PDF of the local hidden state That is

[0191]

[0192] Where the local estimate And the local covariance Are

[0193]

[0194] And, the local prior information for And the local measurement information of neighbors for Can be obtained in the following way

[0195]

[0196] Where

[0197]

[0198] C. Diffusion fusion based on covariance intersection (CI)

[0199] After obtaining the local estimate And the local covariance Of the sensor node s ∈ V, the sensor node s ∈ V exchanges the local estimate and the local covariance with its neighbor nodes, and uses the convex combination linear weighting method to fuse the local estimates and local covariances of the sensor node s ∈ V and its neighbor nodes to update the posterior PDF of the local hidden state That is

[0200]

[0201] Where

[0202]

[0203] Where the weight Can be calculated as follows

[0204]

[0205] The mathematical expectations required for the above variable and parameter updates can be calculated by the following formula:

[0206]

[0207] Step 3: Model Probability Update

[0208] Similar to calculating the model probability in the centralized scenario, in the distributed algorithm, based on Bayes' criterion and the law of total probability, the local model probability of sensor node s can be updated as:

[0209]

[0210] where

[0211]

[0212] where

[0213]

[0214] It can be seen from equations (100) - (105) that the local model probability of sensor node s is updated depending on the corrected local likelihood These corrected local likelihoods can be exchanged among sensor nodes, and then the locally likelihood after diffusion fusion is obtained by taking the convex combination linear weighting method That is

[0215]

[0216] where satisfies and when

[0217] Taking the exponential operation on both sides of equation (106), the locally likelihood after diffusion fusion is

[0218]

[0219] In equation (100), using the locally likelihood after diffusion fusion to replace the original locally corrected likelihood then the local model probability can be calculated as follows:

[0220]

[0221] Step 4: Estimation Fusion

[0222] Similar to the centralized algorithm, the local estimate of sensor node s and can be calculated as

[0223]

[0224] The beneficial effects of the present invention are as follows. The proposed VB-CAIMMF has higher estimation accuracy than the adaptive IMM filter (VB-IMM-QR), the existing interacting multiple model adaptive Kalman filter with a nominal noise covariance matrix (IMM-AKF), and the interacting multiple model intermittent Kalman filter with a true measurement loss sequence and noise covariance matrix (IMM-IKF) under non-ideal conditions such as uncertain, unknown, and time-varying system models, covariance matrices, and measurement loss probabilities. Moreover, in a distributed fusion scenario, the VB-DAIMMF has higher estimation accuracy than Isolated, HCICM, and Diffusion, approaching the optimal estimation VB-CAIMMF. Among them, the Isolated algorithm means that each sensor only uses its own measurement data for state and parameter updates, and there is no information interaction between sensors. The HCICM algorithm means that in the VB-DAIMMF algorithm, a hybrid consistency strategy is selected to replace the CI-based diffusion strategy used in the VB-DAIMMF algorithm for fusion estimation. In this simulation experiment, the number of consistency iterations is selected as 1. The Diffusion algorithm means that a traditional diffusion strategy is used to replace the CI-based diffusion strategy used in the VB-DAIMMF algorithm for fusion estimation. Description of the Drawings

[0225] Figure 1 Schematic diagram of the binary decision factor sequence in a Monte Carlo simulation of Simulation 1;

[0226] Figure 2 Schematic diagram of the RMSE of the position in Simulation 1;

[0227] Figure 3 Schematic diagram of the RMSE of the velocity in Simulation 1;

[0228] Figure 4 Schematic diagram of the true and estimated measurement reception probabilities in Simulation 1;

[0229] Figure 5 ARMSE under different iteration numbers in Simulation 1 pos and running time schematic diagram;

[0230] Figure 6 ARMSE under different iteration numbers in Simulation 1 vel and running time schematic diagram;

[0231] Figure 7 ARMSE under different numbers L of hybrid elements in Simulation 1 pos and running time schematic diagram;

[0232] Figure 8 ARMSE under different numbers L of mixed elements for Simulation 1 vel and the running time schematic diagram;

[0233] Figure 9 Schematic diagram of the sensor network topology for Simulation 2;

[0234] Figure 10 Schematic diagram of the RMSE of the position for Simulation 2;

[0235] Figure 11 Schematic diagram of the RMSE of the speed for Simulation 2. Specific implementation manner

[0236] The present invention will be further described below in conjunction with the accompanying drawings. It should be noted that the following embodiments are based on the present technical solution and give detailed implementation manners and specific operation processes, but the protection scope of the present invention is not limited to this embodiment.

[0237] Embodiment

[0238] The effectiveness and superiority of the proposed VB-CAIMMF and VB-DAIMMF are verified through two target tracking scenarios of single sensor and multi-sensor. The motion trajectory of the maneuvering target can be generated in the following way, that is

[0239]

[0240] wherein, (x 1,k ,x 2,k ) and respectively represent the position and speed of the maneuvering target, Δt = 1s is the sampling interval, Ω is the turning angular velocity, which can be selected from the set {0° / s, 5° / s, -5° / s}, and different turning angular velocities correspond to different motion models. When Ω → 0° / s, this model degenerates into a conventional uniform motion model, that is

[0241]

[0242] The true process noise covariance matrix Q k is set to

[0243]

[0244] where T = 800s, q = 10m 2 / s 3 .

[0245] The model transition probability p ij = 0.90 (i, j ∈ M) and p ij= 0.05 (i ≠ j). The initial model probability is μ0 = [1 / 3, 1 / 3, 1 / 3]. The motion of the maneuvering target is shown in the following table:

[0246] Moment [1s, 100s] [101s, 200s] [201s, 300s] [301s, 400s] Turning angular velocity (° / s) 0 -5 0 5 Moment [401s, 500s] [501s, 600s] [601s, 800s] Turning angular velocity (° / s) 0 -5 0

[0247] The measurement equation where measurement information loss may occur can be expressed as:

[0248]

[0249] where s = 1, 2,..., N s .

[0250] The RMSE and ARMSE of position and velocity are used to evaluate the estimation performance of different algorithms. The RMSE and ARMSE of position are defined as follows

[0251]

[0252] where and are the true and estimated positions of sensor node s at the ι-th Monte Carlo simulation at time k. Mc = 100 is the number of Monte Carlo simulations for the algorithm to run. Similar to the definitions of RMSE pos and ARMSE pos The RMSE vel and ARMSE vel can be calculated by a similar method. To more clearly show the performance of different algorithms, the moving average method with a span of 10 s is used to smooth the RMSE and ARMSE.

[0253] Simulation 1: Single-sensor scenario

[0254] In the single-sensor scenario (N s = 1), the true measurement noise covariance matrix is set to

[0255]

[0256] where r s = 100 m 2 .

[0257] The initial state and the initial covariance P 0|0 are set to

[0258]

[0259] The true and time-varying measurement reception probability of the sensor node is set to

[0260]

[0261] In this scenario, by comparing the proposed VB-CAIMMF(N s = 1) with the adaptive IMM filter (VB-IMM-QR), the interacting multiple model adaptive Kalman filter with a nominal noise covariance matrix (IMM adaptive KF, IMM-AKF), and the interacting multiple model intermittent Kalman filter with true measurement loss sequences and noise covariance matrices (IMM intermittent KF, IMM-IKF). Since IMM-IKF has true measurement loss sequences and noise covariance matrices, in this simulation experiment, it is regarded as a benchmark algorithm for evaluating the performance of other algorithms. In addition, all the above interacting multiple model algorithms use the KLA fusion mode in model condition initialization and estimation fusion. In VB-IMM-QR, the initial nominal state noise covariance matrix and measurement noise covariance matrix are set to and The adjustment parameter τ p = τ r = 500 and the forgetting factor ρ = 1 - exp(-5). In IMM-AKF, the nominal state noise covariance matrix and measurement noise covariance matrix are set to Q k = I4 and R k = 100I2, the initial shape parameter Initial measurement acceptance probability The forgetting factor ρ = 1 - exp(-5). For fairness, in the proposed VB-CAIMMF, the adjustment parameters Initial shape parameter Initial measurement acceptance probability The forgetting factor ρ = 1 - exp(-5), the nominal initial measurement noise covariance matrix The first nominal initial state noise covariance matrix at the initial time The nominal state noise covariance mixing term is set to Initial mixing vector Initial concentration parameter In VB-IMM-QR, IMM-AKF, and VB-CAIMMF, the variational iteration times are set to be the same, all N.

[0262] To better describe the random loss behavior of measurement information in the simulation experiment, Figure 1 The indication factor at all times within 1 - 800 s in a single Monte Carlo simulation experiment for sensor node s = 1 is given Figure 1It can be seen that within 1 - 300 s it indicates that the sensor node s can always accurately receive the measurement data containing useful state information. While in the period of 301 - 800 s, at some moments it indicates that the sensor node s can only collect pure noise information, that is, the useful state information is lost. When N = 5 and the mixing number L = 20, Figures 2 - 3 it shows the RMSE of different algorithms pos and RMSE vel , Figure 4 it gives the true measurement reception probability and the actually estimated measurement reception probability. Table 3 shows the corresponding ARMSE pos , ARMSE vel and the total operation time. From Figures 2 - 3 and Table 4, it can be seen that at the cost of computational efficiency, VB - CAIMMF can obtain better estimation performance than VB - IMM - QR and IMM - AKF, and the estimation performance of VB - CAIMMF is similar to that of the benchmark algorithm IMM - IKF. This indicates that under the conditions of randomly lost measurement data and unknown noise covariance, VB - CAIMMF can still obtain good estimation performance. In addition, from Figures 2 - 3 it can still be seen that within 1 - 300 s, even if there is no loss of measurement data, VB - CAIMMF still has higher estimation accuracy than VB - IMM - QR, which indicates that compared with VB - IMM - QR, the proposed VB - CAIMMF can more effectively reduce the dependence on the initial value of the inaccurate state noise covariance matrix. From Figure 4 it can be seen that the unknown measurement reception probability can be effectively estimated, thus improving the estimation performance of VB - CAIMMF. In the single - sensor scenario, the theoretical computational complexity of VB - CAIMMF can be approximated as while the theoretical computational complexities of IMM - IKF, VB - IMM - QR and IMM - AKF can be approximated as and From the theoretical computational complexities of each algorithm, it can be seen that the computational complexity of VB - CAIMMF is greater than that of other algorithms, which can also be verified in the following table. In addition, it can also be seen that the size of the theoretical computational complexity of VB - CAIMMF depends on {M, N, L}.

[0263] Table 3

[0264] When L = 20 remains unchanged, Figures 5 - 6Shows the ARMSEs of the position and velocity of VB - CAIMMF under different VB iterations and the corresponding operation times. The results show that the estimation performance of VB - CRIMMF gradually improves with the increase of VB iteration times and remains stable when N≥5. When N = 5 remains unchanged, Figures 7 - 8 Shows the ARMSEs of the position and velocity of VB - CAIMMF and the corresponding operation times when L∈[5, 30]. From Figures 7 - 8 It can be found that when L≥5, with the increase of the number of hybrid elements L, the estimation accuracy of VB - CAIMMF will gradually improve, and when L≥20, the estimation accuracy of VB - CAIMMF gradually stabilizes. In addition, from Figures 5 - 8 It can also be found that with the increase of the VB iteration times N and the number of hybrid elements L, the operation time of VB - CAIMMF shows a linear growth trend. Therefore, it is necessary to select appropriate N and L in the actual scenario to balance the relationship between the algorithm estimation performance and the computational complexity.

[0265] Simulation 2: Multi - sensor scenario

[0266] In the multi - sensor scenario, 6 sensors are used to track a two - dimensional maneuvering target, that is, N s = 6. The sensor network topology is as Figure 9 shown. The diffusion weights are determined according to the Metropolis criterion under model j.

[0267] In this scenario, the true measurement noise covariance matrix is set to

[0268]

[0269] where r = 100m 2 .

[0270] For the sensor node s (s = 1, 2,..., 6), its initial state and initial covariance are set to

[0271]

[0272] The true measurement reception probability of the sensor node s (s = 1, 2,..., 6) is set to

[0273]

[0274] In this scenario, the VB-DAIMMF proposed in this chapter is compared with VB-CAIMMF, Isolated, HCICM, and Diffusion. The Isolated algorithm means that each sensor only uses its own measurement data for state and parameter updates, and there is no information interaction between sensors. The HCICM algorithm means that in the VB-DAIMMF algorithm, a hybrid consensus strategy is selected to replace the CI-based diffusion strategy adopted in the VB-DAIMMF algorithm for fusion estimation. In this simulation experiment, the number of consensus iterations is selected as 1. The Diffusion algorithm means that a traditional diffusion strategy is adopted to replace the CI-based diffusion strategy adopted in the VB-DAIMMF algorithm for fusion estimation. In the VB-CAIMMF algorithm, since the data collected by all sensor nodes need to be transmitted to the fusion center for state update and related parameter update, compared with other algorithms (VB-DAIMMF, Isolated, HCICM, and Diffusion), VB-CAIMMF has the optimal estimation performance. Therefore, in this simulation experiment, it is regarded as a benchmark algorithm to measure the advantages and disadvantages of other algorithms. In the above algorithms, all parameter settings are the same, that is: tuning parameter Initial shape parameter Initial measurement reception probability Forgetting factor ρ = 1 - exp(-5), nominal initial measurement noise covariance matrix The first nominal initial state noise covariance matrix at the initial time The nominal state noise covariance mixing term is set to Initial mixing vector Initial concentration parameter The variational iteration number is set to N = 5.

[0275] Figures 10 - 11 Shows the RMSE of different algorithms pos and RMSE vel , the following table shows the corresponding ARMSE pos and ARMSE vel . From Figures 10 - 11As can be seen from the following table, compared with Isolated, VB-DAIMMF, HCICM, and Diffusion, VB-CAIMMF has the best estimation performance, which is reasonable. At the same time, it can also be seen that compared with the Isolated algorithm, the estimation accuracy of VB-DAIMMF in the position and velocity directions has increased by 80.095% and 43.634% respectively. This indicates that the mutual exchange of local information between sensors can greatly improve the estimation accuracy of the algorithm. In addition, it is also found that the estimation performance of VB-DAIMMF is better than that of the Diffusion algorithm, which shows that the mutual exchange of local covariance matrices between sensor nodes helps to further improve the fusion estimation performance. In addition, the single-step running time of VB-DAIMMF is 0.0335 s, while the single-step running time of the Diffusion algorithm is 0.0324 s, which indicates that VB-DAIMMF does not increase more computational burden than the Diffusion algorithm. Compared with the HCICM algorithm, the proposed VB-DAIMMF also has higher estimation accuracy, which may be because the fusion weights used in VB-DAIMMF can be adaptively adjusted according to the local covariance of each sensor node, while the HCICM algorithm uses fixed fusion weights. In addition, the single-step running time of the HCICM algorithm is 0.0313 s, which is lower than that of VB-DAIMMF. This is mainly due to the fact that VB-DAIMMF needs to adaptively calculate the fusion weights.

[0276] Table 4

[0277] For those skilled in the art, various corresponding changes and deformations can be given according to the above technical solutions and concepts, and all these changes and deformations should be included in the protection scope of the claims of the present invention.

Claims

1. An estimation method for distributed multi-model adaptive fusion under unknown noise covariance and measurement loss probability, characterized in that, The estimation method includes: Considering the random loss behavior of sensor measurements, a binary judgment factor obeying Bernoulli distribution is introduced into the sensor node measurement equation, and then the PDF of the sensor node measurement likelihood can be expressed in the form of an exponential product; Within the KLA-IMM framework, by assigning conjugate priors to the unknown noise covariance and the binary judgment factor, the VB technique is used to jointly update the state, noise covariance, and the conditional posterior PDF of the binary judgment factor model; Then, in the mixing and output stages, the KLA method is adopted to fuse the model conditions for estimation, and a VB-based centralized adaptive interacting multiple model filter is obtained; Furthermore, by using a diffusion strategy based on the covariance intersection technique, the obtained VB-based centralized adaptive interacting multiple model filter is extended and applied to the distributed scenario, and a VB-based distributed adaptive interacting multiple model filter is obtained.

2. The estimation method for distributed multi-model adaptive fusion under unknown noise covariance and measurement loss probability according to claim 1, characterized in that, Obtaining a VB-based centralized adaptive interacting multiple model filter includes the following steps: Step 1: Construct a hierarchical Gaussian state space model, that is: Obtain the state transition PDF corresponding to model j, that is Obtain the following exponential product form: Since the binary decision factor follows a Bernoulli distribution, its probability mass function can be written in the following form: Among them, is an unknown and time-varying mixing coefficient, which represents the probability that the measurement data received by the sensor node s at time k has no information loss. Considering that the use of conjugate distributions within the VB framework can greatly simplify the analysis, the beta distribution is chosen to describe the unknown and time-varying mixing coefficients That is Among them, and are shape parameters, For the unknown and time-varying process noise covariance matrix Q k , the unknown state covariance matrix is modeled as a mixture of inverse Wishart distributions as follows: where \(L\) is the number of the mixture of inverse Wishart distributions, and are the degrees of freedom parameter and the scale matrix of the \(l\)-th inverse Wishart distribution corresponding to model \(j\), respectively, \(\chi\) k = [\(\chi\) k,1 ,..., \(\chi\) k,l ,..., \(\chi\) k,L is the mixture coefficient vector, satisfying Introduce an L - dimensional binary random variable where satisfies Then equation (8) can be expressed in the following hierarchical product form Among them, For an unknown and time-varying measurement noise covariance matrix it is modeled using an inverse Wishart distribution, i.e.: Among them, and represent the degree-of-freedom parameter and the scale matrix of the sensor node s, respectively. Considering the slow time-varying characteristics of the mixing coefficient and the noise covariance matrix, a heuristic method can be used to and describe them as follows: where ρ ∈ (0, 1] is the forgetting factor, is the adjustment parameter, the l-th nominal state noise covariance matrix, and is the unknown process noise covariance matrix Q corresponding to model j at time k - 1 k , the mixing coefficient χ k and the unknown measurement noise covariance matrix of sensor node s and the mixing coefficient are the relevant parameters of the posterior PDF Let x k and x k-1 be global latent variables, be the set of common model latent variables and parameters for all sensor nodes in the sensor network, be the set of local latent variables and parameters for each sensor node. The IMM filter obtains the conditional posterior PDF p(Ξ k |m k =j, Z 1:k ) by running a set of sub-filters in parallel corresponding to different models, and obtains the final posterior PDF p(Ξ k |m k =j, Z 1:k ) by further fusing the conditional posterior PDFs p(Ξ k |Z 1:k ). Among them, According to Equations (4)-(11), the joint PDF p(Ξ k , Z 1:k |m k =j) corresponding to model j can be expressed in the following form: wherein, and are the initial state mean vector and the corresponding error covariance matrix corresponding to model j at time k-1, Step 2: Model condition initialization Assume that the posterior PDF corresponding to model i at time k-1 has been obtained, that is According to Bayes' rule and the law of total probability, the mixed probability can be calculated as follows: Similarly, the mixed PDF p (x k-1 |m k = j, Z 1:k-1 ) can be calculated as follows: Similarly, a hybrid PDF can be obtained p (χ k-1 |m k = j, Z 1:k-1 ), p(Q k-1 |m k = j, Z 1:k-1 ) and are as follows: Using the KLA fusion mode, the sum terms on the right-hand side of equations (15)-(19) can be approximated as a single term, that is Among them and can be calculated in the following way: Step 3: Model condition filtering Due to the mutual coupling among the various variables in Ξ k , there is no analytical solution for the conditional posterior PDF p (Ξ k |m k = j, Z 1:k ). To solve this problem, the VB technique is adopted. By choosing a distribution Q j (Ξ k ) corresponding to model j to approximate the conditional posterior PDF p (Ξ k |m k = j, Z 1:k ), based on the mean-field assumption, Q j (Ξ k ) can be expressed as: Use the coordinate ascent method to obtain the optimal solution of the variational distribution q j (Δ), that is Among them At the (t+1)-th (t = 0, 1,..., N-1) iteration, the common model variables and parameters The local variables and parameters of a single sensor node s And the update process of the common hidden state {x k , x k-1} is as follows: A. Common model variables and parameters Update: 1) Let Δ = Q k , substitute Equation (12) into Equation (27) for calculation, then we have where where 2) Let Substitute formula (12) into formula (27) for calculation, then we have where 3) Let Δ = χ k , substitute Equation (12) into Equation (27) for calculation, then we have where B. Local variables and parameters of a single sensor node s (s = 1, 2,..., N s ) Update: 1) Let Substitute Equation (12) into Equation (27) for calculation, then we have wherein, is a normalization constant, and can be obtained by the following formula 2) Let Substitute Equation (12) into Equation (27) for calculation, then we have where 3) Let Substitute Equation (12) into Equation (27) for calculation, then we have where C. Common hidden state {x k , x k-1} Update: 1) Let Δ = x k-1 , substitute equation (12) into equation (27) for calculation, then we have where where 2) Let Δ = x k , substitute equation (12) into equation (27) for calculation, then we have where Among them, and are calculated as follows where In the (t + 1)-th (t = 0, 1, ..., N - 1) iteration process, the required mathematical expectation can be obtained in the following way, and can be calculated as follows: Considering q j,(t+1) (Q k ) and is an inverse Wishart distribution, it can be known that and are Wishart distributions, that is Then and can be calculated as In addition, according to Equations (31) - (38), the required mathematical expectations E j,(t+1) [χ k , E j,(t+1) [logχ k,l , and can be calculated as Step 4: Model probability update Based on the Bayesian criterion, the model probability corresponding to model j can be updated as: where where Step 5: Estimation fusion By using the KLA method to fuse the posterior PDF of the model update state vector, i.e., then and P k∣k can be calculated as Step 6: Initial value selection To implement the proposed VB-CAIMMF algorithm, the initial values of the parameters need to be selected in advance. Assume that the initial mixing coefficient corresponding to model i obeys the beta distribution, that is As can be seen from Equation (64), appropriate shape parameters need to be selected and so as to obtain the prior information corresponding to the initial mixing coefficient of model i For model i, the shape parameters and can be obtained through the mathematical expectation, that is: Among them, is the probability of accurately receiving the initial nominal measurement. Suppose the initial measurement noise covariance matrix corresponding to model i obeys the inverse Wishart distribution, i.e., According to Equation (66), for model i, obeys the Wishart distribution, that is, For model i, the degrees of freedom parameter and the scale matrix can be obtained through the mathematical expectation, that is Among them, is the initial nominal measurement noise covariance matrix of the sensor node s, Similar to Equations (66)-(68), the degrees of freedom parameter and the scale matrix corresponding to the initial process noise covariance matrix Q0 of model i can be obtained through the mathematical expectation of (Q0) -1 That is​ Among them, is the nominal process noise covariance matrix of the first item at the initial moment, Suppose the initial mixture coefficient vector χ0 = [χ 0,1 ,..., χ 0,l ,..., χ 0,L corresponding to model i follows a Dirichlet distribution, i.e., Similarly, for model i, the initial prior concentration parameter can be obtained from the mathematical expectation of the initial mixing coefficient vector χ0, i.e., wherein, is expressed as an initial nominal mixing coefficient vector, and is the l-th initial nominal process noise covariance matrix credibility.

3. The estimation method for distributed multi-model adaptive fusion under unknown noise covariance and measurement loss probability according to claim 2, characterized in that Extending the centralized adaptive interacting multiple model filter to the distributed scenario includes the following steps: Since there is no fusion center in the distributed scenario, local latent variables and local model parameters of each sensor node will be used to replace the global latent variables and common model parameters Similar to the centralized adaptive interacting multiple model filter, local latent variables and local model parameters can be modeled in the following form: Step 1: Model condition initialization In the distributed scenario, a single sensor node s∈V can obtain the conditional posterior PDFs of different models in the same way as VB-CAIMMF, that is Step 2: Model condition filtering In the model conditional filtering stage of the distributed algorithm, VB technology is used to jointly estimate the local hidden state, local variables and parameters. For model j ∈ M, a single sensor node s In the (t + 1)-th (t = 0, 1,..., N - 1) iteration process, it mainly includes three parts, namely: the local hidden state Update of local variables and parameters, local hidden state Update and CI-based diffusion fusion estimation. A. Local state Update of local variables and parameters For sensor node s, its local hidden state Local variables and parameters Posterior Can be updated in a similar way to equations (28) - (43), i.e.: Local process noise covariance Posterior PDF: For model j ∈ M, can be updated to where where Local L-dimensional binary random variable Posterior PDF: For model j ∈ M, can be updated to where Local L-dimensional mixing coefficient Posterior PDF: For model j ∈ M, can be updated to where Local binary decision factor Posterior PDF: For model j ∈ M, can be updated to wherein, is a normalization constant, and can be obtained by the following formula Local mixing coefficient Posterior PDF: For model j ∈ M, can be updated to where Local measurement noise covariance The posterior PDF of: For model j ∈ M, Can be updated to where Local variable Posterior PDF: For model j ∈ M, can be updated to where where B. Local hidden state Update At the (t + 1)-th iteration, a single sensor node s updates the posterior PDF of the local hidden state by using its own local prior information and the local measurement information of its neighbors, i.e., ​ where the local estimate value and the local covariance are Moreover, the local prior information for and the local measurement information for the neighbors can be obtained by the following method where C. Diffusion fusion based on covariance intersection (CI) After obtaining the local estimate of the sensor node s ∈ V and the local covariance the sensor node s ∈ V exchanges the local estimate and the local covariance with its neighbor nodes, and fuses the local estimates and local covariances of the sensor node s ∈ V and its neighbor nodes by means of convex combination linear weighting to update the posterior PDF of the local hidden state i.e., where Among them, the weight can be calculated as follows The mathematical expectations required for the above variable and parameter updates can be calculated by the following formula: Step 3: Model probability update Similar to calculating the model probability in the centralized scenario, in the distributed algorithm, based on the Bayesian criterion and the total probability theorem, the local model probability of sensor node s can be updated as: where where As can be seen from Equations (100)-(105), the update of the local model probability of the sensor node s depends on the corrected local likelihood These corrected local likelihoods can be exchanged between sensor nodes, and then the local likelihood after diffusion fusion is obtained by taking the convex combination linear weighting method That is Among them Meet And when Taking the exponential operation on both sides of Equation (106), the local likelihood after diffusion fusion is In Equation (100), the locally likelihood after diffusion fusion is used to replace the original locally corrected likelihood Then the local model probability can be calculated as follows: Step 4: Estimation fusion Similar to the centralized algorithm, the local estimate of sensor node s and can be calculated as

Images