Comparison of the Extended Kalman Filter and the Unscented Kalman Filter for Magnetocardiography activation time imaging

The non-invasive and radiation-free imaging of the electrical activity of the heart with Electrocardiography (ECG) or Magnetocardiography (MCG) can be helpful for physicians for instance in the localization of the origin of cardiac arrhythmia. In this paper we compare two Kalman Filter algorithms for the solution of a nonlinear state-space model and for the subsequent imaging of the activation/depolarization times of the heart muscle: the Extended Kalman Filter (EKF) and the Unscented Kalman Filter (UKF). The algorithms are compared for simulations of a (6× 6) magnetometer array, a torso model with piecewise homogeneous conductivities, 946 current dipoles located in a small part of the heart (apex), and several noise levels. It is found that for all tested noise levels the convergence of the activation times is faster for the UKF.


Introduction
The localization of the origin of heart arrhythmia is an important part of a successful treatment.For instance, in the context of the Wolff-Parkinson-White syndrome (Nenonen et al., 1991) a pathological accessory pathway is located parallel to the atrioventricular node and can cause serious heart arrhythmia like tachycardia.A non-invasive radiation-free localization of the accessory pathway with ECG or MCG is helpful because it shortens the invasive and X-ray-based catheter mapping procedure.One of the main problems in imaging the electrical activity of the heart is the non-uniqueness of the inverse problem (Fokas et al., 2004) caused by the fact that the number of current dipoles in the heart to be estimated is typically much larger than the number of ECG/MCG sensors.Consequently, the system of linear equations to be solved is ill-posed in general.While in potential imaging the state vector of the state-space model is the transmembrane potential at all heart voxels (Schulze et al., 2009), in activation time imaging the state vector includes the activation times of all heart voxels.The idea behind activation time imaging (He et al., 2002) is to include physiological action potential information, e.g. the wavefront velocity and the upstroke velocity of the depolarization wavefront, to reduce the number of heart model parameters to be estimated without loosing too much accuracy in the calculated sensor signal.The physiological information is incorporated in a cellular automaton model (Weixue et al., 1993) approximating the depolarization wavefront with 3 states and neglecting repolarization.The corresponding nonlinear state-space model can be solved using Kalman Filter algorithms (Liu et al., 2011).It has been shown that activation time imaging on a 3-D myocardium is more stable with respect to measurement noise than potential imaging (Cheng et al., 2003) and can provide an averaged localization error of ∼ 3 mm (Liu et al., 2011) which is small compared to ∼ 1 cm for single dipole localization (Nenonen et al., 1991).In this paper the performances of two different Kalman Filters for the solution of the nonlinear state-space model of activation time imaging are compared for Magnetocardiography: The Extended Kalman Filter and the Unscented Kalman Filter.The convergence of activation times is compared for a quadratic plane (6 × 6) magnetometer array and several measurement noise levels.In Sect. 2 the MCG signal as a function of the activation times is derived, in Sect. 3 the state-space model is explained, and Sect. 4 summarizes the used Kalman Filter algorithms.Section 5 shows the used sensor array, torso model and sources, in Sect.6.1 parameter tests for the Kalman Filters are described, while in Sect.6.2 the EKF and UKF are compared for several measurement noise levels.

From activation times on 3-D myocardium to sensor signal
Following the bidomain theory (Miller and Geselowitz, 1978) we first have to solve the electrical volume conduction problem: with the transmembrane potential V m of the heart muscle cells, primary current density j p , volume current density j V and intracellular conductivity σ i .The volume conduction problem (1) can be solved using the Finite Element Method (FEM) on a 3-D tetrahedral mesh.For the FEM simulation the SimBio neurofem program (SimBio Development Group, 2013) was used.The magnetic field has to be calculated from the total current density j = j p + j V using the Biot-Savart law.For n primary current dipoles J p (t) = (j p (r 1 , t), . . ., j p (r n , t)) T (where the superscript T denotes the transpose) and q magnetometers at a given polarization the superposition principle yields the MCG signal vector caused by the primary currents at the time t: L is the corresponding (q × 3n) leadfield matrix calculated from the volume conduction problem.The cellular automaton approximation for the transmembrane potential at position r in myocardium is shown in Figs. 1 and 2. Now inserting Eq. ( 5) into (2), then Eq. ( 2) into (3) and Eq. ( 3) into (4), and define the activation time state vector as AP approximation without repol.

Fig. 2.
Diffusion of activation in the cellular automaton model (Weixue et al., 1993).τ = 3 ms is the difference of activation times between adjacent voxels.
as well as the corresponding data time series t 1 , . . ., t N .Then the (q × N ) MCG signal matrix φ B = (ϕ B (t 1 ), . . ., ϕ B (t N )) finally is described by where h (a(t − τ (r)) is representing a nonlinear function of τ (r).

State-space model
To tackle the discussed non-uniqueness of Eq. ( 7) the following state-space model is formulated: where x k is the activation time state vector of iteration step k, f is a process function that predicts the k-th state from the (k − 1)-th state, φ B,k is the MCG signal matrix of the k-th heart beat and v k and w k−1 are the measurement noise and process noise, respectively.Both are assumed as Gaussian white noise with the following normal distributions: where R and Q are the noise covariance matrices of the measurement and process, respectively.The process function f is determined by the following rules: where Z is the number of next neighbours (i n ; n = 1, . . ., Z) within a distance d, and v the isotropic activation wavefront velocity.Consequently, the activation time of a voxel tends to the average activation time of the adjacent voxels.In our calculation we chose for the distance to next neighbours d = √ 3d where d is the lattice constant of a cubic grid (Fig. 2).

Kalman Filters
Generally, the application of Kalman Filters is limited to linear processes, and different methods have been proposed to also treat nonlinear ones -as stated by Eqs. ( 9) and ( 7) -by means of this powerful technique.Two of these methods will be presented in the following.

Extended Kalman Filter (EKF)
The idea of the Extended Kalman Filter is the linearization of the nonlinear state-space model Eq. ( 9) by calculating the Jacobian matrices H, F of the functions h, f and approximating the first partial derivatives with difference quotients (Liu et al., 2011): with the temporal resolution ρ and omitting the iteration step index k.The Extended Kalman Filter now minimizes the error covariance matrix defined by: where x true is the true activation sequence and x k the estimated activation of step k.In the "predict" procedure of the algorithm a priori estimates of the state vector and the error covariance matrix are projected from the last step: In the "correct" procedure, the Kalman gain K k is calculated and the state vector and error covariance matrix are updated with the measured data φ B,k : ) x k and P k are then used for the next iteration step.

Unscented Kalman Filter (UKF)
The Unscented Kalman Filter (LaViola, 2003) generates a deterministic set of sampling points, stored in the n×(2n+1) sigma point matrix X k−1 .The columns of X k−1 are calculated by: where ( (n + λ)P k−1 ) i is the i-th column of the matrix square root and λ is defined by: where α and κ are scaling parameters that determine the spread of the sigma points.The square root A of a matrix B satisfies B = AA T and for the symmetric and positive definite matrix B = (n + λ)P k−1 it can be calculated using a Cholesky decomposition (Rhudy et al., 2011).In the "predict" procedure the sigma points are propagated by the process function: Then the a priori state estimate is calculated by: where W (m) i are weights defined by and the a priori error covariance matrix is calculated by with the process error covariance matrix Q and the following weights: β is another scaling parameter to adjust the speed of convergence.In the "correct" procedure first the sigma points are transformed by the measurement function: With the field vector z − k we can now compute the a posteriori state estimate: where the Kalman gain K k of the UKF is defined by with where R is the measurement noise covariance matrix.In the last step the error covariance matrix has to be updated: The computation times of the EKF and the UKF are identical.

MCG sensor array and dipole sources
Figure 3 shows the (6 × 6) array of circular magnetometers with radius r = 9 mm and the torso, lung and heart surfaces used in the simulations.For test purposes the activity is restricted to a small part of the apex of the heart with 946 dipoles distributed on a cubic grid with a lattice constant of d = 1.5 mm.The conductivities of the body tissues were set as in (Liu et al., 2011): torso (0.20 S m −1 ), lungs (0.08 S m −1 ) and cardiac tissue (average of conductivity parallel and transverse to muscle fibre directions: 0.6 S m −1 ).
The non-uniform tetrahedral FEM grid was separated in 48 780 (cardiac tissue), 8045 (lungs) and 238 305 (torso) tetrahedra.The activation wavefront velocity was assumed to be isotropic and we chose the average v = 0.5 m s −1 of the velocity parallel and transverse to the muscle fibre orientation (He et al., 2002) so in the cellular automaton model the difference between the activation times of adjacent voxels is τ = d/v = 3 ms.

Parameter selection, initial state and true state
For the optimization of the convergence of the Kalman Filters several parameters need to be optimized: the scaling parameters α, β, κ for the UKF and the temporal resolution ρ for the EKF.To measure the speed of convergence we define the Relative Error (RE) between the actual state vector x k of step k and the true state vector x true : • 1, process noise covariance Q = 0 and the initial value of the error covariance matrix is the unity matrix (P = 1).Figure 4 shows the parameter test for the temporal resolution ρ of the EKF.The temporal resolution ρ = 3ms with the fastest convergence of the Relative Error is selected for all further EKF calculations.Figure 5 shows the test for the parameter α of the UKF.The other parameters of the UKF are fixed (β = 0 = κ).To enable a stable and fast convergence of the activation times the spread parameter α = 5 was chosen for all further UKF calculations.β and κ have a smaller influence on RE than α.

Comparison of the EKF and the UKF for several measurement noise levels
Figure 6 depicts the convergence of RE as a function of the iteration step k for the EKF (diamonds) and the UKF (dots) for several measurement noise levels with the noise covariance matrix R = σ 2 n •1.The calculated standard deviations of the Gaussian white noise are σ n = 0.01, 0.03, 0.05, 0.1, 0.5 while the minimum of the calculated signal is φ B,min = −3.085and the maximum of the calculated signal is φ B,max = 1.077.For all calculated noise levels Fig. 6 shows that the UKF enables a faster convergence of the RE than the EKF, especially during the first 100 steps.Both for the • 1).Every 5th iteration step is shown.Fig. 6 yields a faster convergence of the activation times of the UKF for all measurement noise levels.The signal and process functions h and f are nonlinear and even non-smooth functions of the activation times.The EKF is an algorithm of 1st order accuracy while the UKF captures the mean and covariance of the Gaussian distributed state vector to the 3rd order Taylor series (Wan and van der Merwe, 2000).In the next steps the heart model needs to be extended to the total ventricles and a muscle fibre model needs to be introduced to include an anisotropic activation wavefront velocity with respect to muscle fibre orientations.For an application to realistic data the heart model parameters, e.g. the process noise covariance Q, need to be estimated from the data.This can be done e.g. with the Expectation-Maximization (EM) algorithm (Khan and Dutt, 2007).

Fig. 1 .
Fig. 1.Approximation of action potential in cellular automaton model: 3 states Resting (R), Wavefront (W) and Excited (E).τ (r) is the activation time at position r in myocardium.

Fig. 3 .
Fig. 3.Left: small heart model with 946 dipoles distibuted on a 3-D cubic grid with lattice constant d = 1.5 mm, Right: (6 × 6) magnetometers with radius r = 9 mm in a quadratic plane array 1 cm above the torso.The magnetic flux is calculated with 16 integration points per sensor.The data sampling rate is t = t j +1 − t j = 3 ms.

Fig. 4 .
Fig. 4. Relative Error RE = ||x k −x true || ||x true || as a function of the iteration step k for the EKF and several temporal resolutions ρ for the difference quotients of the Jacobians H, F (Q = 0, R = 10 −4• 1).Every 5th iteration step is shown.

Fig. 6 .
Fig. 6.Relative Error RE = ||x k −x true || ||x true || as a function of the iteration step k for the EKF (diamonds) and UKF (dots) and several measurement noise levels with noise covariance matrix R = σ 2 n • 1 (Q = 0).Every 5th iteration step is shown.