An Augmented Multiple Imputation Particle Filter for River State Estimation With Missing Observation.

In this article, a new form of data assimilation (DA) method namely multiple imputation particle filter with smooth variable structure filter (MIPF-SVSF) is proposed for river state estimation. This method is introduced to perform estimation during missing observation by presenting new sets of data. The contribution of this work is to overcome the missing observation, and at the same time improve the estimation performance. The convergence analysis of the MIPF-SVF is discussed and shows that the method depends on the number of particles and imputations. However, the number of particles and imputations is influenced by the error difference in the likelihood function. By bounding the error, the ability of the method can be improved and the number of particles and computational time are reduced. The comparison between the proposed method with EKF during complete data and multiple imputation particle filter shows the effectiveness of the MIPF-SVSF. The percentage improvement of the proposed method compared to MIPF in terms of root mean square error is between 12 and 13.5%, standard deviation is between 14 and 15%, mean absolute error is between 2 and 7%, and the computational error is reduced between 73 and 90% of the length of time required to perform the estimation process.

Introduction

The prediction of the river state is important in the hydrology, water resource management, and ecosystem rehabilitation. The knowledge of river flow characteristics is useful for flood forecasting, reservoir operations, and watershed modeling (Sichangi et al., 2018). In flood forecasting, the predicted state is use to produce alerts of the incoming flood to prevent damages to human life, properties and environment (Jain et al., 2018). Besides that, the information from the state estimation is used in controlling the outflows of reservoir during low flows day of river and also rapid flows resulting from dam-break that may cause catastrophy to the environment and massive losses to life and property (Adnan et al., 2017), (Cao et al., 2019). The water flow can also be used in watershed modeling that is used to manage the water by channeling the water from any sources into a single larger body of water such as a larger river (Aswathy et al., 2016). The management and planning of the water source is important for the increasing of the water demand in the next few years due to the increasing population growth, urbanization, industrial use, irrigation needs, and water-intensive agriculture (Tinka et al., 2013). This includes the planning of the water projects, irrigation systems, hydropower system, and optimized utilization of water resources (Adnan et al., 2017). The information of the river discharge (flow) is necessary in the management procedure and the variation of the hydrologic cycle is related to the climate change, land use, and water use (Bjerklie et al., 2018). The river discharge is one of the climate variables by the Global Climate Observing System (Tarpanelli et al., 2019). The prediction of the variables that represent the water flow regime helps in the ecosystem rehabilitation program that seeks to safeguard or restore indigenous ecosystems by manipulating the river flow (Blythe and Schmidt, 2018). The estimation of the river state can be carried out using the data assimilation (DA) method. The DA method is a mathematical technique that combines observation data with the system model and creates the updated model state while maintaining the parameter of the model (Li et al., 2014). The updated model state is defined as the probability density function (pdf) that is based on Bayes’ theorem and known as the new posterior pdf (Smith et al., 2013). The new state is acquired whenever there are new observations and uses them to initiate the next model forecast as reported in a study (He et al., 2014). The DA method is desired to optimally and consistently estimated, even if the noisy readings arrive sequentially in time (Liu and Gupta, 2007). The DA technique is accessible in two classes namely, variational and sequential (Gadsden et al., 2014).

The variational method is based on the optimum theory of control. Optimization is carried out on the associated parameters by minimizing the cost function that influences the model to misfit the information. The examples of this technique are variational data assimilation (VAR), evolutionary data assimilation (EDA), and maximum likelihood ensemble filter (MLEF) (Abaza et al., 2014). The VAR method makes estimates by minimizing the cost function which measures the difference between the model estimate, the observation, and the associated uncertainties. The gradient-based optimization algorithm is used to adjust the model states and parameters of the model and receive the appropriate estimate for the measurement (Kim et al., 2014). Besides that, the EDA utilizes a multi-objective evolutionary strategy to continually develop the set of model states and parameters where the model error and penalty function minimization for each assimilation time step is determined adaptively. This is done to enhance the convergence of parameters and lead to excellent estimation outcomes (Bertino et al., 2003; Solonen et al., 2014). Another type of variational method is the MLEF, which combines the VAR and the ensemble Kalman filter (EnKF). This method maintains the strength of both method and capable to handle nonlinear model dynamics as well as nonlinear observation. However, in some cases the performance of the MLEF may deteriorate due to consistency of the observation (Rafieeinasab et al., 2014). Besides that, the sequential methods use a probabilistic framework and estimate the whole system state sequentially by propagating information only forward in time. This method does not require an adjoint model and makes it easy to adapt with the model (Rigatos, 2012). The sequential-based technique is frequently used in estimation compared to the variational method, since the prior state is updated with the new observation available and the process is performed sequentially (Arulampalam et al., 2002). The examples of this method are extended Kalman filter in (Li et al., 2014), ensemble Kalman filter (EnKF) in (Rafieeinasab et al., 2014), unscented Kalman filter (UKF) (Pintelon et al., 2017), cubature Kalman filter (Liu and Gupta, 2007), and particle filter (PF) (Ugryumova et al., 2015). The EKF is commonly used sequentially DA method due to easy implementation, but depends strongly on the accuracy of the system linearization that is performed using Taylor series expansion (Li et al., 2014). This technique shows excellent efficiency in low nonlinearity and diverges in higher nonlinear cases (Zhang et al., 2015). For a highly nonlinear system, the estimation can be carried out using the EnKF that offers estimation without linearization (Rafieeinasab et al., 2014). This method involves a large sample size to represent the number of samples or ensemble members for precise estimation that can be produced with Monte Carlo method, Latin hypercube sampling, and moment equation (Kang, 2013; Gadsden et al., 2014; Wang et al., 2017). The samples’ mean and covariance are used to perform updating in the estimation process (Rafieeinasab et al., 2014). However, the large sample size may cause for computational demand (Gadsden et al., 2014). Another form of EnKF is the ensemble square root filter (EnSRF) that does not require the observation to be perturbed as the standard implementation of the EnKF. The algorithm of this method has been demonstrated to be as fast as the EnKF and more precise for a specified ensemble size than EnKF (Liu et al., 2017).

The UKF is more appropriate for application for low computational and high accuracy. This method includes the sigma points derived from the unscented transformations (Ding et al., 2015). The sigma points are propagated through the system model and the related weight factor during the estimation process and produce new sets of sigma points, which are subsequently used in the calculation of the projected states (Pintelon et al., 2017). However, this method relies on the precise of prior noise distribution. Wrong prior value can lead to large estimation errors or divergence of errors (Mao et al., 2017). Another type of sequential DA method is the PF, which does not involve system linearization and includes a number of particles during prediction (Hu et al., 2012). Each particle represents the estimated state with its associated likelihood that is determined using the residual between the simulated output and observation (Solonen et al., 2014). Large numbers of particles provide the chance of good estimation, but the computational time will increase. Another factor to consider when using PF is the issue of degeneration and dimensionality curse that can influence the estimation (Ugryumova et al., 2015). All the abovementioned DA methods were implemented in the river state estimation study. The modification of the methods may improve the performance of the method and the result obtained. In a study (Cao et al., 2019), the particle filter is modified by considering the system condition and the particle weighting procedure. The modification has improved the estimation result. Besides that, the modification can also be made by combining two or more methods such as the combination of the PF with smooth variable structure filter (SVSF). The SVSF produce estimated state and state error covariance, which is used to formulate the proposal distribution for the PF to generate particles (Ogundijo et al., 2016). The SVSF is robust to modeling errors and uncertainties (Feng et al., 2011). This method produces better estimated result than the PF only. Besides that, the SVSF also combined with the cubature Kalman filter (CKF). The CKF offers the nearest known approximation to the Bayesian filter in the sense of maintaining second-order information contained in the Gaussian assumption of noisy measurements. This method also does not need Jacobians and therefore applies to a wide range of issues. The accuracy of CKF and the stability of the SVSF ensure good estimation results (Liu and Gupta, 2007).

In this work, the estimated river system flow and stage, and also the velocity of the sensor are inspired from (Tinka et al., 2013) and (Zhang et al., 2014). These values can be obtained using the previously mentioned DA method. However, the missing observation data can be problematic that will affect the estimation process (Ismail and Jalaludin, 2016). The missing data can be handled using modern method or traditional method. The modern method is represented by the maximum likelihood and multiple imputations, while the traditional techniques are deletion and mean imputation techniques (Gadsden et al., 2012). It is noted that the multiple imputation particle filter (MIPF) is introduced to deal with this problem by randomly drawn values known as imputations to replace the missing information and then uses the particle filter to predict nonlinear state as reported in (Habibi, 2007). However, the addition of the new data burdens the estimation process. Therefore, the SVSF is introduced to bound the error difference and assist in the estimation process.

In this article, the convergence analysis of the MIPF–SVSF in estimating the river flow and stage from the initial condition is presented. The article is constructed as follows: the system model, observation model, and the state space model for estimation process are presented in Modelling the Marine River Flow. Problem Formulation briefly explains the effect of the missing observation during estimation. Then, the proposed algorithm for estimation with missing data is described in the Proposed Data Assimilation Approach. The convergence analysis is presented in Convergence Analysis that includes almost sure convergence and convergence of the mean square error. In Results and Discussion, the details on estimation process and numerical simulations are discussed and finally conclusions are presented in Conclusion.

Modeling the Marine River Flow

The river flow model can be represented by one or two-dimensional Saint–Venant equations (Tinka et al., 2013) depending on the characteristics of the water flow. If the flow is in one-dimensional, the 1-D Saint–Venant equations is considered. However, if the flow is not one-dimensional, which may happen in flood plains or in large rivers, the 2-D Saint–Venants equation is more suitable to be applied (Litrico and Fromion, 2009). Besides that, the representation of the observation is referring to the movement of the sensor since the Lagrangian sensor is used in this research (Tinka et al., 2013). The combination of the system model and the observation is represented by the state-space model and used in the DA method.

System Model and Observation Model

By considering one-dimensional flow of the river without any uncontrolled release of water flow, the system model is represented by 1-D Saint–Venant equations. These two equations coupled are first order hyperbolic partial differential equations (pde) derived from the conservation of mass and momentum. By considering a prismatic channel that have same cross-section throughout the length of channel with no lateral inflow, the equation is represented as (Tinka et al., 2013). where A is the cross section, Q is the discharge or flow, L is the river reach, T is the free surface width, D is the hydraulic depth, S is the friction slope, S is the bed slope, F is the gravitational acceleration, h is the distance of the centroid of the cross section from the free surface, P is the wetted perimeter, and m is the Manning roughness coefficient.

Observation Model

The system observation is represented by the velocity of the flow as measured by the sensors. The relation between the velocity of the sensor and the flow at the corresponding cross-section relies on assumptions made about the profile of the water velocity that considered as the observation model. The profile is the combination of the average velocity in the transverse and vertical direction. In transverse direction, the surface velocity profile is assumed to be quartic, and the von Karman logarithmic profile is assumed in the vertical direction. By considering a particle moving at a distance y from the center line and z from the surface, the relation between the particle’s velocity and the water flow is represented by the following equations (Tinka et al., 2013) with where w is the channel width; d is the water depth; A , B , and C are constants; and K is the Von Karman log constant.

State Space Representation

During estimation process, the system and observation equation is represented as a state-space model that comprises of the parameters of the model, observation, system noise, and measurement noise. The channel is discretized into n cells with each cell and has same length. In the nonlinear system state estimation, the initial conditions and the boundary conditions of the system are required as the inputs. The uncertainties of the model and also the inaccuracies of the inputs measurements are considered as the system noise . While the measurement noise, represent the errors and uncertainties of the measurements. Both noises are represented by the zero mean Gaussian error. The state-space model for the estimation is described as follows where X is the state vector at time t and the input u contains the boundary conditions, i.e., the upstream flow and downstream stage. where and are the flow and stage at cell i at time t, respectively, and n is number of cells used for the discretization of the channel. Since the system is observed by n sensors, Eq. 10 can be reformulated into where denotes the noisy observation of the state X such that the is an independent and identically distributed (i.i.d.) measurement noise and g is the measurement transformation for sensor k.

Problem Formulation

The Bayesian theorem used by the DA method is represented as follows (Crisan and Doucet, 2002). where is the system state at time t, is the observation at time t, is the posterior probability of state, at time t given observation, from time one to time t, is the likelihood function of state at time t given observation at time t, is the prior probability of state at time t given observation from time one to time t-1, is the normalizing constant.

In this theorem, the observation data is used to adjust the likelihood function. The adjusted likelihood function is used to modify the prior probability to obtain the desired posterior probability that represents the estimated state. The normalizing constant in this theorem is represented as (Ugryumova et al., 2015). where all parameters are defined in Eq. 14. Based on Eqs 14, 15, the posterior probability is very much depending on the likelihood function . This function is represented as follow (Mukherjee and Sengupta, 2010) where and are defined in Eq. 14, and represent the estimation of the observation using the predicted states. The difference between the observation and the estimated observation is considered in this function. This error affects the likelihood function, and thus influences the prediction result.

In this research, several types of observation namely y and z positions of the sensors, and the velocity of the sensors are considered. The positions of the sensors are used in determining the estimated observation using Eqs 4–8. Next, the obtained estimated velocity is compared with the measured velocity of the sensors and form the likelihood function for this case. Since the likelihood function is important in estimation, the continuous observations from the sensors are desired to secure this function throughout the estimation process. In the event of missing observation data, the likelihood function is affected, and thus limits the ability of the standard DA method. Therefore, the MIPF method is introduced to perform estimation with new input data that replace the missing data. The availability of the observations is checked at each time instance. The missing data are handled by introducing a random indicator variable, R (Abaza et al., 2014). Consider the overall observation comprises both available and missing data from all sensors. The observation at time instance t for all sensors k = 1…,n with is defined as the missing information set , while the available information set is the observation for all sensors k = 1,…, n such that .

The introduction of the new data may affect the error difference between the observation from the new data and the estimated observation. Therefore, the SVSF method that is robust and stable in the estimation process is introduced to handle this problem. The combination of the MIPF and SVSF capable of handling state estimation with missing information and error differences problem.

During missing information, several random observations or imputations is introduced in the estimation process. The imputations are drawn from the proposal function (Chai and Draxler, 2014) where represent all missing observations at time t, represent all available observation from time 0 to time t, is the particle set with no regard of missing data, N is the total number of particles and M is the total number of imputations. Next, the imputations are reformulated into the imputed data sets of where represent all missing observation during jth imputation and time t, and represents all available observation at time t. The imputed sets are used in determining the posterior probability density that represent as where is the system state at time t, is the complete observation, and are defined in Eq. 17, and is defined in Eq. 18. Considering the Monte Carlo approximation, the probability density can be written as where is defined in Eq. 17, is defined in Eq. 18, and , , and , are defined in Eq. 19. For each data set , the probability density from particle filtering is written as follows: where is the system state at ith particle and jth imputation at time instance t, is the related weight. By substituting Eq. 21 into Eq. 20, the overall representation of the desired posterior probability density is represented as where , M, N are defined in Eq. 17, is defined in Eq. 19, and , are defined in Eq. 21. A smooth variable structure filter (SVSF) is a sliding mode-based predictor-corrector estimator. By having the proper representation of the switching gain, the estimation is converged to be within the boundary of the true state values. This ensures that the estimator is stable and robust to modeling uncertainties and noise. The width of the boundary is referring to the existence subspace that represents the number of uncertainties present in the estimation process. The uncertainties are associated with the inaccuracy of the internal model of the filter and measurement model, that is varies with time. The selection of the width is based on a prior knowledge, since the parameter is not exactly known. During estimation process with proper representation of the boundary, the estimated states are forced to switch back and forth along the true state trajectory by the SVSF gain. However, the switching may cause for chattering effect and can be reduced by introducing the smoothing subspace. For smoothing subspace that is bigger than the existence subspace, the chattering effect is reduced. While for smaller smoothing subspace compared to the existence subspace, the chattering effect is still present (Feng et al., 2011).

Proposed Data Assimilation Approach

The algorithm of the MIPF–SVSF is the combination of the MIPF and SVSF. The MIPF is functioning to handle the nonlinearity and missing data problem, while the SVSF is used to deal the noise problem. The algorithm of the MIPF–SVSF method is presented as follows:

Initialization

where is the initial state, is the initial prediction error covariance, is the system noise covariance matrix, is the measurement noise covariance matrix, and is the initial weight for particle filtering and is the initial error difference between the observation and the estimated measurement.

Prediction

Draw random imputation/observation from proposal function (as in MIPF) from the previously available measurement where is the missing observation at time t, represent all available observation from time 0 to time t, is the system state during available observation, is the weight for particle filtering during available observation, is the overall measurement that includes the available observation and new imputation, N is the number of particles, and M is the number of imputations. where is the input, is the previous system state, is the one step estimated state, is the estimated measurement, is the previous prediction error covariance, is the new prediction error covariance, is the system noise covariance, and A is the linearized system model by using Jacobian as represented in Eq. 27. where is the convergence rate, is the error difference between the overall observation (that includes the new observation that replaces the missing observation and the available observation) and the estimated measurement, is the linearized measurement model by using Jacobian as represented in Eq. 31, is the smoothing boundary layer vector, and is the diagonal matrix constructed from number of represented as the saturation function, is defined by:

- Prediction in SVSF involves producing the one step state estimation and the covariance P for particles generation

- Updating in SVSF involve updating the state estimation and covariance P for particles generation

- Generate N particles of M imputation by using the previously obtained estimated state, and prediction error covariance,

Updating

Determine the important weight of the particles from the previous weight, and the likelihood function that includes the generated particles, and the overall measurement . The important weight, is normalized with the sum of particles weight is equal to unity Next, the effective sample size, from the weights are calculated, in order to measure the degeneracy problem. If the is smaller than the threshold (e.g.: <0.5xN), severe degeneracy problem might occur. So, the particles with small weights are eliminated and concentrate on the particles with large weights. The eliminated particles are replaced with the new set of particles from resampling process and the weight is represented as in Eq. 41. where is defined in Eq. 38, is the new weight if the degeneracy problem occur, N is defined in Eq. 14. By using the generated particles and their associated weights, the estimated states are represented as follows where is defined in Eqs 39, 41, is defined in Eq. 37. For the purpose of next iteration in estimation, the updated measurement estimate, is determined and used to produce the updated measurement error, as shown in Eqs 43, 44 respectively. where is defined in Eq. 25, is defined in Eq. 42.

Convergence Analysis

In order to perform the convergence analysis, the state-space model of the system and observation, and the MIPF–SVSF are reformulated into probability representation.

Probability Space Formulation

Let be a probability space where is the Borel set of , the Borel set is the standard set of all possible probability events on . Two types of vector-valued stochastic process namely system state, and observation, are involved in this space. The system state, X is a Markov process of initial distribution X ∼µ and probability transition kernel, K(x |x ). where , , n is the dimension of the states, n is the number of the states. Since there are three states namely flow, stage, and cross section, the current states are represented as , and the states at previous time, . The overall observation is represented as for of n sensors that is independent to each other and have marginal distribution where is the overall observation by kth sensor at time t, is the dimension of the observation and is defined in Eq. 45. For missing observation problem, consider the non-response vector-valued stochastic process, with dimensional vector. The availability of the observation is indicated using and introduced the following sets where is the missing observation at time t, is the available observation at time t, is the indicator for the availability of the observation, and are defined in Eq. 46. The probability density of the non-response mechanism that is corresponding to the proposal function to draw imputation as in Eq. 17 is represented as where is the missing observation at time t, is the available observation from time 0 to time t, and is the indicator for the availability of the observation from time 0 to time t, and are defined in (17).

Probability Representation for State Estimation

The estimation using MIPF–SVSF is represented by the posterior probability density function that considers both the available observation and the missing observation. The distribution of the overall posterior probability density and the probability density of the states during available observation, are described as follows: where represent the available observation from time to time , is the time for state during the available observation, is defined in Eq. 17, and are defined in Eq. 45, , , and are defined in Eq. 48. For notational convenience, the probability density is written as . So, the distribution of and are related as follows: where is the probability density of the states at time t, and are defined in Eq. 48, is defined in Eq. 49. By applying the standard Bayesian filtering theory, the posterior probability density of the available observation can be expressed as where is the probability density of the states at time t-1, is defined in Eq. 51, is defined in Eq. 45, is defined in Eq. 46.

By combining Eq. 51 with Eq. 52, the overall posterior probability density can be reformulated into Eq. 53. This equation shows the relationship between the probability density during available observation, and the overall probability density, . Since the missing observation is very much affecting the probability distribution, the additional knowledge of this problem is covered by introducing the empirical distribution, to (51) as a replacement to the true distribution of the probability density, . So, the posterior probability density for each state with N particles is represented as where is the available observation at time t, is the probability density of each states during available observation, is defined in Eq. 48. The empirical distribution, consists of N particles that are distributed approximately according to . given is the set of particles distributed approximately according to for all state, given by as in the standard bootstrap procedure. Next, the missing observation is incorporated into Eq. 32 by applying naive Monte Carlo approximation with for . The posterior probability density of the desired states with N and M is represented as where is the posterior probability density for the related states with N particles, is the missing observation at time t and jth imputation, N and M are defined in Eq. 17. By referring to the MIPF–SVSF algorithm, the relation between the approximated distributions of the related states and the true density can be analyzed through convergence analysis.

Part 1: Prediction Using SVSF

Consider the one step prediction error where is the estimate of at time with the initial state, is the one step prediction at time , and are defined in Eq. 1.

By using the Taylor series expansion around , is linearized as follows where , and is the high order terms of the Taylor series expansion where is a problem dependent scaling matrix, is introduce to provide extra degree of freedom to the filter, is an unknown time-varying matrix for linearization errors of the dynamical model that satisfies that gives where is defined in Eq. 59, is defined in Eq. 1, is defined in Eq. 60, and are defined in Eq. 60. Applying the Taylor series expansion around to the , the innovation of the filter is represented as where , is a problem dependent scaling matrix, is introduce to provide extra degree of freedom to the filter, is an unknown time-varying matrix for linearization errors of the dynamical model that satisfies as in Eq. 61, is defined in Eq. 62, is defined in Eq. 59, and are defined in Eq. 2. The filtering error is represented as where is the gain, is the estimated state from SVSF, is defined in Eq. 45, is defined in Eq. 59, are defined in Eq. 63, is defined in Eq. 2,

The one step prediction error covariance, is given by where is the previous prediction covariance, is defined in Eq. 23, is defined in Eq. 59, is defined in Eq. 1, is defined in Eq. 60, and are defined in Eq. 60

The filtering error covariance is given by where is defined in Eq. 23, is defined in Eq. 64, are defined in Eq. 63.

Given matrices A, H, E, and F with appropriate dimensions such that . Let X be a symmetric positive definite matrix and be an arbitrary positive constant such that . Then the following inequality holds

Theorem 3: Consider Theorem 1 and Theorem 2 and assume and as in Eqs 60, 62 are true. Let be positive scalars. The upper bound for the one step prediction error covariance matrix and filtering error covariance matrix can be represented by the following Riccati-like difference equations using Lemma 1: where is the gain like in Eq. 64, and are defined in Eq. 23, is defined in Eq. 59, is defined in Eq. 60, is defined in Eq. 60, are defined in Eq. 63, are positive scalars.

With initial condition have positive definite solutions and such that for all the following two constraints are satisfied.

The filter gain for the upper bound is referring to the in Eq. 32 and minimizes the upper bound so that

For , suppose that and . If and Then the solutions for to the following difference equations are Satisfy.

By referring to Lemma 2, rewrite the error covariance matrices as the function of . where is defined in Eq. 65, , are defined in Eq. 23, is defined in Eq. 59, is defined in Eq. 1), , , and are defined in Eq. 60, is defined in Eq. 56, are defined in Eq. 63.

Next, rewrite the upper bound obtained in Eq. 67 and Eq. 68 as the function of and where is defined in Eq. 61, and are defined in Eq. 23, is defined in Eq. 59, , is defined in Eq. 60, are defined in Eq. 63, are positive scalars.

Consider Eqs 75–78, the error covariance matrices and the upper bound satisfy condition Eq. 72, 73 in Lemma 2. which gives .

Part 2: Prediction Using MIPF

Consider: if is a sequence of random probability measures, then converges to if for any continuous bounded function

The likelihood function is a bounded function in argument and represented as .

For any with random variables obtained from states, and prediction error covariance matrices, from SVSF as in Eq. 37. Then, for any where is the probability density of the random variables, is the probability density of the one step state estimate, and are the constant, is the particles obtained, and is the number of the states; N is described in Eq. 37.

Consider Let be the generated by , then and, as Using Minkowski’s inequality where .

Part 3: Updating Using MIPF

For any

Then for any where is the probability density of the one step state estimate, is the probability density of the state without consideration of missing data, and are the constant, N is described in Eq. 37.

Proof: Consider For the first part of Eq. 90 Using Minkowski’s inequality

Part 4: Updating With New Imputation Using MIPF

For any where is the probability density of the state without the consideration of missing data, is the probability density of the state with the consideration of missing data, and and are the constant; M and N are described in Eq. 37.

. Then, for any that includes new imputation

Proof: Consider For the first part of Eq. 95 For the second part of Eq. 94 that refer to Eq. 89 Using Minkowski’s inequality

Part 5: Resampling Using MIPF

For any and There exists a constant , such that for any and where is the probability density of the state with the consideration of missing data, is the posterior probability density of the states, and and are the constant; M and N are described in Eq. 37.

Consider the following By using Minkowski’s inequality Let be the generated by , then and That gives

Theorem 4For all there exist that is independent of N but being influenced by M for any bounded function, where is the posterior probability density of the states, is the constant, M and N are described in Eq. 37, is particles generated as in Eq. 37. The number of imputations, M for missing observation data is related to the number of particles, N since at each imputation contain several particles for estimation. To ensure the accuracy of estimation, the number of imputation M and the number of particles N are considered based on the dimension of the missing observation data, and dimension of the states, respectively. Besides that, the number of M and N also are influenced by the Gaussian noise of the measurement error that affects the likelihood function as in Eq. 16. Large error difference in the likelihood function requires high number of M and N to help for the convergence of the filter to the true state. By having proposed states and error covariance matrices to generate the particles as in Eq. 37 that are determined through considering the error difference of the observation and estimated measurement, the distribution of the particles can be limit to the relevant area only. Therefore, the number of M and N needed by this method will be less than the particle filtering method that randomly generated particles without reference.

Results and Discussion

In this section, the estimation of the river flow and stage, and the velocity of the last drifter are presented. Consider the measurements are suffering from the missing data that affect the estimation process. The MIPF–SVSF is proposed to perform the state estimation by applying several numbers of particles and imputations as described in the previous section. The performance of this method is evaluated by finding the root mean square error (RMSE), standard deviation (SD) (Chai and Draxler, 2014), mean absolute error (MAE), and computational time (Liang et al., 2012) between the measured velocity and the estimated velocity. The RMSE, SD, MAE, and computational time of the MIPF–SVSF are compared with the same parameters from forward simulation, EKF during complete data, PF during complete data, and MIPF.

Description of the Estimation Process

In this research, the measurements from the drifters (Lee et al., 2011) are used in the estimation of the flow, stage, and cross section that are later used to estimate the velocity of the sixth sensor. The estimated sensor velocity is proportional to the velocity of the river flow. During estimation process, the river system is discretized into 60 cells with 5 m interval. The temporal step size is chosen as 1 s. Since the data about the bottom of the river is unavailable, the bed slope is considered as zero.

The estimation is carried out by combining the system and the observation. During no missing data, the estimation is carried out by the PF–SVSF. The MIPF–SVSF is proposed for the state estimation with missing data problem. Consider the observation is suffering from three types of missing data namely missing velocity data, missing position data, and missing combination of velocity and position data at the same time. Whereby each case has 10, 20, and 30% missing data (Ugryumova et al., 2015). During estimation, the missing data is replaced by several number of imputations. The imputations are generated based on the previously available data. For each imputation, the error difference between the estimated measurement by the SVSF, and the observation from the imputation is calculated. This error is bounded by using the smoothing boundary layer vector. The bounded error and the convergence rate are used in producing the SVSF gain to obtain the new state value.

Next, the new state by the SVSF with their covariance are used to generate several particles. The particles and the observation from the imputation is combined to form the likelihood function. From the likelihood function, the important weight of the particles is obtained. The particles and their weight are used to produce estimated state that consists of the flow, stage, and cross-sectional area. The number of estimated states is depending on the number of imputations. The mean of these states is the desired new estimated state by the proposed method. The estimation process is repeated for 400 s with the error difference between the current result and the observation is use in the next estimation cycle. Since the proposed method includes the convergence rate and smoothing boundary layer vector, the error difference between the estimated velocity and the observation is reduced. By reducing the error, the MIPF–SVSF will have less computational time compared to the MIPF.

Discussion

The estimation process by the DA method includes the merging of the system and the observation. Three types of missing data cases are considered in this research as mentioned earlier. For each case, several numbers of imputations, namely 5, 10, 15, and 20 are injected during estimation with 50 particles. The same number of particles are applied for the whole estimation process due to the good estimation result by the PF using these particles during no missing data. Figure 1 shows the estimated flow and stage by forward simulation, the EKF (no missing data), MIPF and MIPF–SVSF for 30% missing velocity data. The overall flow estimation by the MIPF is smaller than the EKF and forward simulation, while the overall estimation by the MIPF–SVSF is smaller than the other methods. Besides that, the MIPF also produce estimated stage that is bigger than the EKF and the forward simulation. However, the estimated stage by the MIPF–SVSF is bigger than the MIPF. The state estimation by the proposed method seems to be reasonable. Therefore, the performance of the method is evaluated by finding the velocity of the final drifter and compared with the measurement. The velocity of the final drifter is obtained by combining the estimated flow and cross-sectional area at the corresponding cell. Figure 2 shows the velocity of the final drifter predicted by the forward simulation, EKF (no missing data), MIPF, and MIPF–SVSF. The figure shows that the proposed method helps the estimation to converge to the desired value. This can be seen by smaller difference between the estimated and the true state, compared to the other methods.

FIGURE 1

Estimated flow and stage at the 30th cell, forward simulation, EKF during no missing data, MIPF with 30% missing data, and MIPF–SVSF with 30% missing data.

FIGURE 2

Estimated velocity of the final drifter using forward simulation, EKF during no missing data, MIPF with 30% missing data, and MIPF-SVSF with 30% missing data.

For further analysis and comparison of the performance of the MIPF–SVSF, the RMSE, the SD, the MAE, and the computational time are determined as in Table 1, Table 2 and Table 3. The performance of the MIPF–SVSF is compared with the MIPF and the PF and the EKF from Table 4. Since the proposed method is the upgraded version of the MIPF, the performance of the MIPF is first examined. By considering all missing data cases, the MIPF is able to produce RMSE, SD, and MAE that is close to the estimation result by the PF during complete data. These show that the external input data successfully fill in the missing part. Besides that, the combination of the SVSF with the MIPF improves the estimation result through smaller RMSE, SD, and MAE compared to the MIPF for all missing data cases. This indicates that the SVSF method successfully bound the error difference in the estimation process. The bounded error also helps to reduce the MIPF–SVSF’s computational time compared to the MIPF through the generation of gain that brought the estimation velocity closer to the observation. The results are represented by bold numbers in Tables 1–3. The selection is made based on the smallest value among variables.

TABLE 1

Performance of the MIPF and MIPF–SVSF during missing velocity data.

Method	N	M	Missing velocity data
			10%				20%				30%
			RMSE	SD	MAE	Time	RMSE	SD	MAE	Time	RMSE	SD	MAE	Time
MIPF	50	5	0.500	0.437	0.344	59.158	0.500	0.439	0.351	60.564	0.507	0.444	0.351	61.945
		10	0.502	0.434	0.345	61.363	0.498	0.437	0.349	63.936	0.508	0.440	0.348	66.452
		15	0.506	0.439	0.348	62.623	0.500	0.441	0.347	67.522	0.507	0.441	0.346	71.077
		20	0.506	0.436	0.351	64.263	0.503	0.438	0.346	70.894	0.502	0.438	0.346	75.787
MIPF-SVSF	50	5	0.434	0.373	0.322	6.809	0.445	0.376	0.338	8.047	0.443	0.371	0.334	8.890
		10	0.437	0.376	0.333	8.233	0.443	0.376	0.333	10.892	0.451	0.379	0.332	13.102
		15	0.435	0.372	0.329	9.716	0.434	0.365	0.327	13.703	0.445	0.379	0.333	16.891
		20	0.435	0.375	0.330	11.221	0.438	0.372	0.329	16.682	0.442	0.373	0.331	20.699

Bold values represents the smallest value among variables (RMSE SD MAE Time).

TABLE 2

The performance of the MIPF and MIPF–SVSF during missing position data.

Method	N	M	Missing position data
			10%				20%				30%
			RMSE	SD	MAE	Time	RMSE	SD	MAE	Time	RMSE	SD	MAE	Time
MIPF	50	5	0.493	0.432	0.347	64.212	0.497	0.432	0.347	66.177	0.501	0.435	0.348	67.914
		10	0.500	0.432	0.346	65.994	0.502	0.434	0.348	70.1.1	0.502	0.433	0.345	73.766
		15	0.503	0.434	0.344	68.312	0.503	0.435	0.344	74.529	0.503	0.435	0.345	79.754
		20	0.502	0.435	0.347	70.403	0.502	0.433	0.348	77.947	0.503	0.434	0.348	84.730
MIPF–SVSF	50	5	0.439	0.362	0.338	7.162	0.441	0.373	0.333	8.494	0.441	0.371	0.336	9.155
		10	0.440	0.369	0.335	8.778	0.438	0.372	0.333	10.972	0.441	0.376	0.333	13.482
		15	0.442	0.372	0.333	10.239	0.442	0.373	0.333	13.979	0.440	0.366	0.336	17.643
		20	0.441	0.382	0.334	11.705	0.440	0.370	0.334	17.022	0.441	0.370	0.336	21.743

Bold values represents the smallest value among variables (RMSE SD MAE Time).

TABLE 3

The performance of the MIPF and MIPF-SVSF during missing combination of position and velocity data.

Method	N	M	Missing combination of velocity and position data
			10%				20%				30%
			RMSE	SD	MAE	Time	RMSE	SD	MAE	Time	RMSE	SD	MAE	Time
MIPF	50	5	0.502	0.436	0.346	65.196	0.504	0.436	0.343	66.356	0.502	0.438	0.351	67.764
		10	0.502	0.435	0.347	67.258	0.502	0.436	0.348	70.214	0.504	0.435	0.346	73.130
		15	0.502	0.436	0.345	68.972	0.501	0.433	0.346	75.441	0.501	0.433	0.348	76.950
		20	0.501	0.437	0.349	71.112	0.497	0.435	0.347	77.984	0.499	0.438	0.352	82.081
MIPF–SVSF	50	5	0.438	0.372	0.330	6.866	0.443	0.374	0.337	7.960	0.446	0.384	0.334	8.594
		10	0.441	0.378	0.332	8.231	0.441	0.374	0.333	10.821	0.442	0.374	0.335	12.233
		15	0.444	0.377	0.337	9.647	0.438	0.371	0.331	13.794	0.448	0.379	0.334	15.838
		20	0.442	0.374	0.336	11.175	0.439	0.376	0.330	16.579	0.436	0.375	0.324	19.397

Bold values represents the smallest value among variables (RMSE SD MAE Time).

TABLE 4

The performance of the Forward simulation, EKF, and PF during complete data.

Method	N	RMSE	SD	MAE	Time
Forward sim	—	1.6778	1.054	1.417	1.511
EKF	—	0.657	0.526	0.480	2.735
PF	50	0.504	0.436	0.339	58.272

By using 50 particles with different number of imputations, the performance of state estimation by MIPF–SVSF during missing velocity data is shown in Table 1. The result shows that the increasing number of imputations will reduce the RMSE, SD, and MAE. However, too many imputations may diverge the estimation from the true value. This can be seen through the sudden increase of RMSE, SD, and MAE after reduction. Besides that, the number of imputations also related to the percentage of missing data. Whereby the higher percentage of missing observation data require higher number of imputations. The same response is shown during the state estimation with missing position data and missing combination of velocity and position data as shown by Tables 2, 3, respectively. The increasing number of imputations reduces the RMSE, SD, and MAE, and the suitable number of imputations is selected before the divergence of the estimation occurred. For different percentages of missing data, the appropriate number of imputations is considered based on the response from RMSE, SD, and MAE.

In this research, the number of particles is fixed to 50 particles throughout the estimation process. Based on the result represented by bold numbers in Tables 1–3, the MIPF–SVSF has less RMSE, SD, MAE, and computational time compared to the MIPF for all missing data cases. The percentage improvement of the proposed method compared to MIPF in terms of RMSE is between 12 and 13.5%, SD is between 14 and 15%, and MAE is between 2 and 7%. While for the computational time, the percentage of reduction is between 73 and 90%. Therefore, the proposed method still can have small RMSE, SD, MAE and computational time under lower number of particles. Since the number of particles and imputations are related, a smaller number of imputations for MIPF–SVSF compared to the MIPF is enough for good estimation performance.

Conclusion

The MIPF–SVSF is the extension of the MIPF method with the addition of the SVSF. The proposed method introduces several numbers of imputation based on the previously available data to replace the missing data. There are three types of missing data considered in the research, namely, missing velocity data, missing position data, and the missing combination of velocity and position. The convergence analysis of this method shows that the number of particles and imputations depends on the likelihood function that represents error difference between the estimated observation and the imputation data that replaces the missing data. Large error difference requires high number of particles and imputation to converge the estimation to the true state. The SVSF reduces the error difference through the introduction of convergence rate and the smoothing boundary layer vector. These variables are used to generate the particles and their weight and form the estimated state. The performance comparison between MIPF–SVSF and MIPF shows that MIPF–SVSF has better performance than the MIPF in terms of RMSE, SD, MAE and computational time. Besides that, different missing data cases with different percentages of missing data require different combination of particles and imputations. The MIPF–SVSF requires smaller combination of particles and imputation compared to the MIPF for all missing data cases.