Conception and analysis of Cascaded Dual Kalman Filters as virtual sensors for mastication activity of stomatognathic craniomandibular system

The presented work shows a possible methodical approach for parameter estimation of a kinematic and dynamic element that characterizes a human mandibular system during the mastication process using position measurement only. The considered parameters are the velocity, friction coefficient, and the mass of the moving part of the mandibular during the mastication activity of a human. Internal or optical motion sensors can still allow imprecision in the measurements. To overcome these, in the present work a system identification algorithm is designed using a combination of three backward cascaded Kalman Filter, which consists of three Extended Kalman Filters. The identification procedure is validated through a matching criterion based on the estimation of the mass, which is assumed to be known in the first stage of the Kalman Filter structure. Three EKFs are tuned as long as the initial value of the mandibular mass is achieved as an estimation of the third one. This is due to the fact that the optimization procedure tries to optimize a non-convex optimization problem that can admit more than one solution. The main contribution of this project is designing state estimation dynamic system, which accurately estimates friction with a linear time varying model. Friction coefficient plays an important role in the early diagnosis of temporomandibular joints disorders, since it is very low under normal condition, and an increase may be associated with abnormalities. Computer simulations show the effectiveness of the proposed method to accurately estimate friction dynamics and refrain from complex nonlinearities.


Introduction
In the fields of virtual sensing, sensor fusion, intelligent sensing and sensorless control, the Kalman Filter (KF), with its methodological variations, is one of the highly known algorithms used for linear systems. The Extended Kalman filter (EKF) can be used for the identification of parameters of nonlinear systems. Today sensor fusion plays an important role in motion analysis in all types of applications. The recording of motion is often done by inertial or optical motion sensors. However, the accuracy of such tools is often very low, which can result in incorrect data interpretation. Kalman Filter is used as a real time multi-sensor fusion algorithm to overcome these problems, see [1]. Additionally, it serves for estimating the states of a system [2], to increase the accuracy [3] or to make corrections for system parameters [1]. Sequential filter such as KF "provides an excellent framework for estimating and updating the parameters of a system as data becomes available" [4], therefore a dual KF structure is proposed. The recording of the lower jaw movements and positions under neuromuscular and kinematic view serves: • for the qualitative and quantitative determination of the functional state of the the position and movement of the jaw; • to gain indications for therapeutic measures (function-oriented therapy planning); • for the documentation of the results due to the course of therapy in order to control and evaluate the initiated therapy (functional therapy documentation and feedback therapy).
Its purpose is the improvement of function in terms of a biomechanical optimisation or prosthetic rehabilitation and/or oral rehabilitation on a function-oriented basis. The data collected from the movement records represent a function-oriented therapeutic approach. They provide patient-oriented values for setting the articulator for the design of individual occlusal surfaces. In the digital workflow, the instrumental motion analysis is an important and necessary supplement in order to map the movement functions of the jaw for the occlusal design of restorations in the CAD/CAM process. Human oral cavity consists of both hard and soft tissues such as chin, teeth, tongue, mucosa and glands as well as the temporomandibular joint (TMJ) which connects the upper and lower jaw bone. The TMJ can be considered as a synovial joint and works in the same way as other synovial joints such as the hip and the knee. The loading and motion conditions in synovial joints are complex therefore "it is important to measure the friction in synovial joints accurately, since the friction level in these natural bearings is generally low and also because of additional difficulties associated with other soft tissues surrounding the joint and mechanical factors that can contribute to the friction measurement" in the context of [5]. The main motivation of this project is to estimate the friction which is generally very low in biological systems under normal conditions, but can drastically increase in deteriorated conditions or due a disease. The main function of the oral cavity as an organ is closely related to speech and food processing. Therefore, friction plays an important role. [6] studies possible effects of high friction and TMJ disc disorders. According to them, for example, due to some changes in the lubrication mechanism the coefficient may rise by a factor of about 4. In [6] it is employed friction as a Coulomb-type friction model, which, however, is considered as a simplification when taking into account characteristics of friction dynamics between a deformable and a rigid articular surfaces. Coulomb-type friction model is inherently nonlinear and relies upon time and velocity dependence of friction [7]. Another well-known dynamical friction estimation model is the LuGre model. In [8], designing a model-based nonlinear controller called Dynamical Adaptive Backstepping-Sliding Mode Control (DAB-SMC) for a human like skeletal-muscle the authors implemented LuGre-based adaptive friction observer. With the latter controller gains more stability in performance, and other metrics, such as position tracking accuracy, are also improved. However, in order LuGre-based friction model, for instance, to incorporate hysteresis behavior of friction in a more complex model structure, the parameter estimation complexity increases due to its high nonlinearity [8]. In this contribution a form of a linear time varying friction model is proposed. In fact, the friction parameter is considered in an EKF and dynamically estimated. This approach is a valid compromise between a structurally simple model and a dynamic estimation of the parameter. This approach can be seen as a linearization of a complex phenomenon such as the mechanical friction effect around the actual working state (position and velocity) of the mandibular system. The effectiveness of the obtained results indicates that the adopted method is a suitable one for the considered application.
"KF is used for estimating an object's position and other parameters of motion based on measurements, which are non linear related to the objects position" [9]. For this reason, the goal was to use the advantages of the KF also for the approximation of a jaw position. The paper presents three EKFs to be tuned as long as the initial value of the mandibular mass is achieved as an estimation of the last EKF. This is due the fact that here in the model it tries to optimize a non-convex optimisation problem which can admit more than one solution. The paper is structured in the following way. The paper starts with the description of the physical model in Section 2. It follows a description of the Kalman Filter background in Section 3 and the Observability of the system in Section 4. The main results are presented in Section 5, with a presentation of the used measurement system and the proposed algorithm to estimate velocity and parameters of the moving mandibular system. A discussion of the found results will be addressed in Section 6. The conclusion closes the paper.

Physical model
The concept of using kinematics in translational and rotational capacity is used for the combination of opening-closing and lateral movements of the jaw in [10]. Therefore, the calculated forces, combined with the moment arms and horizontal components can be used for the calculation of the forces that contribute to the lateral movements. The following equations describing the Kinematic Model of the jaw joint system are proposed in [10]. The opening and closing movement M n can be described as where n = LM, RM, LD or RD: left masseter (LM), right masseter (RM), left diagastric (LD) and right diagastric (RD) describes the product of the force component in the sagittal plane, M A n describes the moment arm and a subject-specific parameter K n is added to account for inaccuracies. Therefore the total joint moment about the centre of rotation M tot can be determined using where W is the mandible weight, D is the joint damping and O z is the offset value that determines resting position in the vertical direction. The kinematics of the joint can be calculated by where d(T s + ∆T s ) describes the new angular displacement (in θ-direction) of the mandible, on the previous displacement d(T s ), previous angular velocity ω(T s ), moment of intertia I and sampling time ∆T s . In order to obtain vertical displacement z(T s ) and not polar coordinates, a geometric conversion that uses (R, θ) was used in [10]. The vertical displacement z(T s ) is given by The transitional component of the movement can be determined similarly, but along the x-axis.
The horizontal force components can be calculated using the acute angle between force vector and x-axis, and is described by α. The vertical distance y can be calculated by α can be determined as follows with F x being the lateral component of the total musculotendon force F The total force F Tot , can be determined by the sum of the forces in x-direction, with the addition of O x to add lateral force offset and subject specific parameters K i (i=LD, RD, KM, RM) to account for different muscle sizes and strengths The horizontal displacement x(T s + ∆T s ) and v(T s + ∆T s ) can be calculated as follows The horizontal motion also requires its own damping components F Damp due to the passive viscous and elastic damping properties of the masticatory system that only affect horizontal movement where β is the damping coefficient and v the horizontal velocity. The adopted model is based on the following mathematical systems: including the three states of the system in which z mp (t) describes the position of the mandible, z mv (t) the velocity of the mandible,ż mv (t) the acceleration of the mandible, and u(t) is total force. The parameter that is being estimated by the system varies. In this specific case the parameter k v , which describes the Friction coefficient is being estimated.We getk v (t) = 0 as it is set as stochastic constant. The Jacobian matrix of the system is as follows Estimating the Mass m the system can be described as follows therefore the Jacobian matrix can be defined as 3. Kalman filter background "The Kalman Filter is a set of mathematical equations that provide an efficient computational (recursive) means to estimate the states of a process, in a way that minimizes the mean of the squared error." [11] Based on [11] the process of a non-linear system can be described with the following stochastic difference equation with a measurement defined as in which w k and v k both represent random variables describing the process and measurement noise respectively. Therefore, the a priori estimation of the state is defined aŝ as well as the a priori error covariance matrix P − k , which is as follows The matrix Q k is the process noise covariance and both J k and W k are the process Jacobians at the current step k. Therefore P k−1 represents the a posteriori error covariance matrix at step k − 1 (previous step). Secondly the EKF measurement update equations include the Kalman gain K which is as follows with H k and V k as the measurement Jacobians and R k as the measurement noise covariance all changing at each time step k with the sampling time T s . The system which is being considered in this paper neglects both W k and V k , since the noises are assumed to have the same effect on all states and measured variables and are therefore not considered further. The a posteriori estimation is as followsx where z k (k) = h(x − k , 0). The a posteriori estimation of the error covariance matrix is as follows In which H is the Jacobian matrix of the output measurement.

Observability of the system
The analysis of the observability is done for each singular system. For that purpose, the idea is to consider separately in a cascaded form a procedure for the estimation of the mastication parameters. Considering the full system with four states implies the following The observability analysis can start with a observable determination matrix θ as follows: only the position is a measured state, then It follows that To have a sufficient condition for the observability, then θ = n. The observability test for system (24) follows as: and This is already proven that the sufficient condition for the full system is not guaranteed. Considering the following system note for the kinematic Forward Euler model to estimate the state of the force F Tot (k − 1) and acceleration γ(k − 1): Calculating as before, then θ 0 = 1 0 0 , and The sufficient condition is guaranteed when Rank θ = 3. Now, considering the system of estimating the Friction coefficient k v represented with the following Kalman filter: where u(t) is already estimated F Tot ,ż mv (t) is the already estimated velocity, which we use together with z mp (t) as measured states: Given two measured states the observable determination matrix (27) can be obtained as: and Hence, above described system is always observable iff z mv (t) = 0 and u(t) = 0. The system for the estimation of the mass m is defined as follows: The same procedure of calculation of the corresponding elements of observability determination matrix follows as: and Based on the rank criteria of observability it is easy to show for θ matrix with above described elements that Rank θ = 3 holds true. Thus, the observability of discussed system (40) is guaranteed.

Experimental setup
The JMA-Optic System (Zebris Medical GmbH, "Isny", Germany), shown in Fig. 1a, which is a 3D coordinate measuring system, was used to measure the movement of the mandible. It consists of three parts: a handy stand-alone face bow which contains a receiver and a control unit, the lower jaw sensor which works as the transmitter unit and a charging station. The system works due to infrared signals from the mandible sensor when measuring from the facebow. During movement cameras in the headbox are detecting a change of the image of the mandible sensor. The coordinates of the mandible sensors are being calculated in 3D space.
A patented bite fork, shown in Fig. 1b is used for an exact relationship between the movement data in the measuring system and the position of the upper jaw as well as the surfaces of the teeth scanned by the model and intraoral scanner. It allows an easy transfer of the settings of the jaw to mechanical articulators. The measured position data consist of three fixed points on the bite fork. One point is located in the coordinate origin, the other two are equally spaced to the right and left of it. The position values of the three points were queried at the same frequency during the motion sequence.
A simplified mathematical model of the jaw joint system was designed. It includes parameters such as a combination of all acting forces, described as force d, the acceleration a of the movement of the mandible, a spring coefficient k s , a friction coefficient k v and the mass m of the mandible. As it has already been described in the previous sections, it is impossible to determine the exact mass of the lower jaw, therefore it is usually assumed that the entire mass is applied to a centered point and only acts at this point. The measured position data was used to determine the parameters mentioned above. For this the motion sequence along the y-axis (vertical axis) was used. A cascaded KF structure was created for this purpose. As there is no input a dual KF was used. It was not possible to approximate all parameters at the same time, as the system would have been unobservable otherwise. Based on the cascaded structure, the parameters could thus be continuously estimated and at the same time compared and checked against given reference values from the literature. Due to the relative large sample time (0.0167 sec.) with which the position is sampled, the dynamics of a feedforward discretisation appears to be high and, thus, makes it difficult to be tuned. Figure 2 represents a systematic construction of the structure and indicates also the proposed estimation procedure. The three EKFs are tuned as long as the initial value of the mandibular mass are achieved as an estimation of the third one. The diagram shows that the cascaded structure consists of three individual EKF's, with different known inputs and different output. Each EKF can therefore determine an individual parameter (EKF 1 estimates velocity z mv (k), whereas EKF 2 estimates k v ). The EKF's are separate from each other, otherwise the system would be unobservable.

Estimation of friction coefficient
Given are the equations of the EKF 2 as already generally described in (13) which approximates k v . After applying Backward Euler discretisation the system structure can be written as the following, where T s represents a sampling time: x − (k) describes the a priori estimation of the three states x 1 (k), x 2 (k) and x 3 (k), in which where H describes the output Jacobian Matrix. The Jacobian matrix of the system is as follows Ts 1+Ts

Estimation of mass
Given are the equations of the EKF 3 as already described in (15) which estimates m.
1+Ts kv(k)/x 3 (k) + in which x 1 (k − 1) is the position, x 2 (k − 1) the velocity and x 3 (k − 1) the mass m. The Jacobian matrix for it can be described as while the output Jacobian matrix follows as:

Results
The determination of the parameters was done by means of the cascaded EKF structure, which allowed a continuous approximation and feedback comparison of the parameters. As shown in Fig. 2, the mass in EKF 1 is assigned a fixed value. In the last step of the cascaded structure (EKF 3 ) the mass is then being estimated as an unknown parameter and compared with the given value. Thus it can be determined whether the estimates of whole cascaded structure are valid. The measured position and the estimated position by EKF are presented in Fig. 3 with a  Figure 4 shows the measured velocity and the estimated velocity by EKF. Only small deviations are visible in the approximation of measured velocity. Figure 5 shows the approximation of the friction coefficient k v . According to reference value [12] estimated k v of temporomandibular joints varies between approximately 0.021 −0.023 Ns/m. As we can see, for the whole simulation friction value demonstrates steady approximation of these values. Figure 6 shows the estimation of the mass m which approaches a constant value of approximately 0.25 kg. The estimation of the mass is used as a validation approach. The mass was known in EKF 1 and EKF 2 (Fig. 2) Figure 4: Result of the measured velocity and estimated velocity by EKF it is emphasized that the other parameters such as k v , total force as well as the velocity which are being estimated during the course of the cascaded structure approached are reliable values, as the estimation of the mass includes the previously estimated other parameters. However, the tuning of the KF must be optimal for the generation of such constant and comparable results. Therefore, in each EKF of the estimation of friction coefficient and mass a tuning of Q x was individually applied in compliance with cascaded structure, and R w was manually calibrated to achieve the best empirical approximation.

Conclusions
The presented topic represents an interesting and important approach for dentists as well as for the industry, which deals with the required measuring systems. Observing friction of temporomandibular joints in the oral cavity can be an important indicator for early disease diagnosis considering crucial functionality of the oral cavity in speech and food processing. Based on demonstrated simulation results the proposed backward cascaded Kalman Filters  Figure 6: Approximation of the mass of the mandible which consist of three backward Extended Kalman Filters is a promising to identification algorithm, where under adverse abnormal condition friction coefficient tends to deviate from inherent biological low level. The results showed that for incomplete observable system designed cascaded structure still allows to accurately estimate through a validation process of target mass. The choice of backward discretisation due to the large sampling time of the position effectively restrained the effect of high dynamics.
Further research work will focus on the improvement of the discussed cascaded Kalman Filters structure in terms of reducing the deviation of real states from filtering results while incorporating more measurements. At the same time, for the simulation an adaptive filtering process should be considered to adjust and update the Kalman gain based on adaptive tuning of noise covariance matrices, which will also presumably lead to an improved EKF robustness.