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).

$T\frac{\partial H}{\partial t}+\frac{\partial Q}{\partial x}=0(1)$

$\frac{\partial Q}{\partial t}+\frac{\partial }{\partial x}\left(\frac{Q^2}{A}\right)+\frac{\partial }{\partial x}\left(g{h}_{c}A\right)=gA\left(S_o-S_f\right)(2)$

$S_f=\frac{m^2{Q}^2{P}^{\frac{4}{3}}}{A^{\frac{4}{3}}}(3)$

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_f is the friction slope, S_o is the bed slope, F_g is the gravitational acceleration, h_c 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)

$v_p\left(y,z\right)=F_T\left(y\right)F_v\left(z\right)\frac{Q}{A}(4)$

with

$F_T\left(y\right)=A_q+B_q{\left(\frac{2y}{w}\right)}^2+C_q{\left(\frac{2y}{w}\right)}^4(5)$

$A_q+B_q+C_q=0(6)$

$A_q+\frac{B_q}{3}+\frac{C_q}{5}=1(7)$

$F_v\left(z\right)=1+\left(\frac{0.1}{K_v}\right)\left(1+\log \left(\frac{z}{d}\right)\right)(8)$

where w is the channel width; d is the water depth; A_q, B_q, and C_q are constants; and K_v 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 $v_t$. While the measurement noise, ${\epsilon }_t$ 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

$X_{t+1}=f\left(X_t,u_t,{\upsilon }_t\right)(9)$

$W_t=g\left(X_t,{\epsilon }_t,t\right)(10)$

where X_t is the state vector at time t

$X_t={\left(Q_2^t,\dots ,Q_{n_c}^t,H_1^t,\dots ,H_{n_c-1}^t\right)}^T(11)$

and the input u_t contains the boundary conditions, i.e., the upstream flow and downstream stage.

$u_t={\left(Q_1^t,H_{n_c}^t\right)}^T(12)$

where $Q_i^t$ and $H_i^t$ are the flow and stage at cell i at time t, respectively, and n_c is number of cells used for the discretization of the channel. Since the system is observed by n_k sensors, Eq. 10 can be reformulated into

$W_t=\left(\begin{array}{c}{g}_1\left(X_t,{\epsilon }_{t,1},t\right)\\ ⋮\\ g_{n_k}\left(X_t,{\epsilon }_{t,n_k},t\right)\end{array}\right)=\left(\begin{array}{c}{W}_{t,1}\\ ⋮\\ W_{t,n_k}\end{array}\right)\text{for}k=1\dots n_k(13)$

where $W_t$ denotes the noisy observation of the state X_t such that the ${\epsilon }_{t,k}$ is an independent and identically distributed (i.i.d.) measurement noise and g_k 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).

$p\left(X_t|W_{1:t}\right)=\frac{p\left(W_t\text{|}{X}_t\right)p\left(X_t\text{|}{W}_{1:t-1}\right)}{p\left(W_t\text{|}{W}_{1:t-1}\right)}(14)$

where $X_t$ is the system state at time t, $W_t$ is the observation at time t, $p\left(X_t|W_{1:t}\right)$ is the posterior probability of state, $X$ at time t given observation, $W$ from time one to time t, $p\left(W_t\text{|}{X}_t\right)$ is the likelihood function of state $X$ at time t given observation $W$ at time t, $p\left(X_t\text{|}{W}_{1:t-1}\right)$ is the prior probability of state $X$ at time t given observation $W$ from time one to time t-1, $p\left(W_t\text{|}{W}_{1:t-1}\right)$ is the normalizing constant.

In this theorem, the observation data $W_t$ 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).

$p\left(W_t|W_{1:t-1}\right)={\displaystyle \int p\left(W_t|X_t\right)p\left(X_t|W_{1:t-1}\right)d{X}_t}(15)$

where all parameters are defined in Eq. 14. Based on Eqs 14, 15, the posterior probability is very much depending on the likelihood function $p\left(W_t\text{|}{X}_t\right)$. This function is represented as follow (Mukherjee and Sengupta, 2010)

$p\left(W_t\text{|}{X}_t\right)=p\left(W_t-g\left(X_t\right)\right)(16)$

where $W_t$ and $X_t$ are defined in Eq. 14, and $g\left(X_t\right)$ 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_t,k (Abaza et al., 2014).

$R_{t,k}\{\begin{array}{l}0:\text{Observation}\text{is}\text{missing}\text{from}\text{sensor}k\text{at}\text{time}t\\ 1:\text{Observation}\text{is}\text{available}\text{from}\text{sensor}k\text{at}\text{time}t\end{array}$

Consider the overall observation $W_{t,k}$ comprises both available and missing data from all sensors. The observation at time instance t for all sensors k = 1…,n_k with $R_{t,k}=0$ is defined as the missing information set ${\Xi }_t$, while the available information set ${\text{Ψ}}_t$ is the observation for all sensors k = 1,…, n_k such that $R_{t,k}=1$.

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)

${\Xi }_t^j\sim \phi \left({\Xi }_t\text{|}{\Psi }_{0:t}\right)={\displaystyle \sum _{i=1}^N{\tilde{\omega }}_t^{i}p\left({\Xi }_t|{\tilde{X}}_t^i\right)}\text{for}i=1,\dots ,N\text{and}j=1,\dots ,M(17)$

where ${\Xi }_t$ represent all missing observations at time t, ${\text{Ψ}}_{0:t}$ represent all available observation from time 0 to time t, ${\left\{{\tilde{\omega }}_t^i,{\tilde{X}}_t^i\right\}}_{i =1}^N$ 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

$U_t^j=\left\{{\Xi }_t^j,{\Psi }_t\right\}(18)$

where ${\Xi }_t^j$ represent all missing observation during jth imputation and time t, and ${\Psi }_t$ represents all available observation at time t. The imputed sets are used in determining the posterior probability density that represent as

$p\left(X_t|{\Psi }_{0:t}\right)={\displaystyle \int p\left(X_t|W_{0:t-1},{\Psi }_t\right)p\left({\Xi }_t|{\Psi }_{0:t}\right)d{\Xi }_t}(19)$

where $X_t$ is the system state at time t, $W_{0:t-1}$ is the complete observation, ${\Xi }_t$ and ${\text{Ψ}}_{0:t}$ are defined in Eq. 17, and ${\text{Ψ}}_t$ is defined in Eq. 18. Considering the Monte Carlo approximation, the probability density can be written as

$p\left(X_t|{\Psi }_{0:t}\right)\approx \frac{1}{M}{\displaystyle \sum _{j=1}^{M}p\left(X_t|W_{0:t-1},U_t^j\right)}(20)$

where $M$ is defined in Eq. 17, $U_t^j$ is defined in Eq. 18, and $X_t$, ${\text{Ψ}}_{0:t}$, and $W_{0:t-1}$, are defined in Eq. 19. For each data set $U_t^j$, the probability density from particle filtering is written as follows:

$p\left(X_t|W_{0:t-1},U_t^j\right)\approx {\displaystyle \sum _{i=1}^N{\omega }_t^{j,i}\delta \left(X_t-X_t^{j,i}\right)}(21)$

where $X_t^{j,i}$ is the system state at ith particle and jth imputation at time instance t, ${\omega }_t^{j,i}$ is the related weight. By substituting Eq. 21 into Eq. 20, the overall representation of the desired posterior probability density is represented as

$p\left(X_t|{\Psi }_{0:t}\right)\approx \frac{1}{M}{\displaystyle \sum _{j=1}^M{\displaystyle \sum _{i=1}^N{\omega }_t^{j,i}\delta \left(X_t-X_t^{j,i}\right)}}(22)$

where ${\Psi }_{0:t}$, M, N are defined in Eq. 17, $X_t$ is defined in Eq. 19, and $X_t^{j,i}$, ${\omega }_t^{j,i}$ 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

${\widehat{x}}_0,Cov{P}_{0|0},CovR,CovQ,{\omega }_0,e_0(23)$

where ${\widehat{x}}_0$ is the initial state, $Cov{P}_{0|0}$ is the initial prediction error covariance, $CovQ$ is the system noise covariance matrix, $CovR$ is the measurement noise covariance matrix, and${\omega }_0$ is the initial weight for particle filtering and $e_0$ 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

${\text{Ξ}}_t^j\sim \varphi \left({\text{Ξ}}_t|{\text{Γ}}_{0:t}\right)={\displaystyle \sum _{i=1}^N{\tilde{\omega }}_t^{i}p\left({\text{Ξ}}_t|{\tilde{x}}_t^i\right)}fori=1,\dots ,Nandj=1,\dots M(24)$

$W_t^j=\left\{{\text{Γ}}_t,{\text{Ξ}}_t^j\right\}(25)$

where ${\Xi }_t$ is the missing observation at time t, ${\Gamma }_{0:t}$ represent all available observation from time 0 to time t, ${\tilde{x}}_t^i$ is the system state during available observation, ${\tilde{\omega }}_t$ is the weight for particle filtering during available observation, $W_t^j$ 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.

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

${\widehat{x}}_{t|t-1}=f\left({\widehat{x}}_{t-1|t-1},u_{t-1|t-1}\right)(26)$

$A_t=\frac{\partial f}{\partial x}{|}_{{\widehat{x}}_{t|t-1}{u}_{t-1|t-1},}(27)$

$Cov{P}_{t|t-1}=A_t.Cov{P}_{t-1|t-1}.A_t^T+Cov{Q}_{t|t-1}(28)$

${\widehat{\zeta }}_{t|t-1}=g\left({\widehat{x}}_{t|t-1}\right)(29)$

where $u_{t-1|t-1}$ is the input, ${\widehat{x}}_{t-1|t-1}$ is the previous system state, ${\widehat{x}}_{t|t-1}$ is the one step estimated state, ${\widehat{\zeta }}_{t|t-1}$ is the estimated measurement, $Cov{P}_{t-1|t-1}$ is the previous prediction error covariance, $Cov{P}_{t|t-1}$ is the new prediction error covariance, $Cov{Q}_{t|t-1}$ is the system noise covariance, and A is the linearized system model by using Jacobian as represented in Eq. 27.

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

$e_{W,t|t-1}^j={\text{W}}_t^j-{\widehat{\zeta }}_{t|t-1}(30)$

$B_t=\frac{\partial g}{\partial x}{|}_{{\widehat{x}}_{t|t-1},}(31)$

$K_{svsf,t}^j=B_t^{-1}.diag\left[\left(\left|e_{\text{W},t|t-1}^j\right|+\vartheta \left|e_{\text{W},t-1|t-1}^j\right|\right).sat\left({\text{Ψ}}^{-1}{e}_{\text{W},t|t-1}^j\right)\right]\cdot diag{\left(e_{\text{W},t|t-1}^j\right)}^{-1}(32)$

${\widehat{x}}_{t|t}^j={\widehat{x}}_{t|t-1}+K_{svsf,t}^j{e}_{\text{W},t|t-1}^j(33)$

$Cov{P}_{t|t}^j=\left(I-K_{svsf,t}^{j}B\right)Cov{P}_{t|t-1}^j\left(I-K_{svsf,t}^{j}B\right)T+K_{svsf,t}^{j}Cov{R}_{t|t}{\left(K_{svsf,t}^j\right)}^T(34)$

where $\vartheta$ is the convergence rate, $e_{\text{W},t|t-1}^j$ 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, $B_t$ is the linearized measurement model by using Jacobian as represented in Eq. 31, $\Psi$ is the smoothing boundary layer vector, and ${\Psi }^{-1}$ is the diagonal matrix constructed from $n_k$ number of $\Psi$ represented as

${\text{Ψ}}^{-1}=\left[\begin{array}{ccc}\frac{1}{{\text{Ψ}}_1}& 0& 0\\ 0& \ddots & 0\\ 0& 0& \frac{1}{{\text{Ψ}}_{n_k}}\end{array}\right](35)$

the saturation function, $sat\left({\Psi }^{-1}{e}_{\text{W},t|t-1}^j\right)$ is defined by:

$sat\left({\text{Ψ}}^{-1}{e}_{\text{W},t|t-1}^j\right)=\left\{\begin{array}{c}1\\ e_{{\text{W}}_{\text{i}},t|t-1}^j/{\text{Ψ}}_{\text{i}}\\ -1\end{array} \begin{array}{c},e_{{\text{W}}_{\text{i}},t|t-1}^j/{\text{Ψ}}_{\text{i}}\ge 1\\ ,-1<e_{{\text{W}}_{\text{i}},t|t-1}^j/{\text{Ψ}}_{\text{i}}<1 \\ e_{{\text{W}}_{\text{i}},t|t-1}^j/{\text{Ψ}}_{\text{i}}\le -1\end{array}\right\}fori=1,\dots ,n_k(36)$

- Generate N particles of M imputation by using the previously obtained estimated state, ${\widehat{x}}_{t|t}^j$ and prediction error covariance, $Cov{P}_{t|t}^j$

$x_t^{j,i}\sim q\left({\widehat{x}}_{t|t}^j,Cov{P}_{t|t}^j\right)fori=1,\dots ,Nandj=1,\dots ,M(37)$

Updating

Determine the important weight of the particles from the previous weight, ${\omega }_{t-1}^{j,i}$ and the likelihood function that includes the generated particles, $x_t^{j,i}$ and the overall measurement $W_t^j$.

${\widehat{\omega }}_t^{j,i}\sim p\left(W_t^j\text{|}{x}_t^{j,i}\right).{\omega }_{t-1}^{j,i}(38)$

The important weight, ${\widehat{\omega }}_t^{j,i}$ is normalized with the sum of particles weight is equal to unity $\left({\displaystyle \sum _{i=1}^N{\widehat{\omega }}_t^{j,i}}\right)$

${\omega }_t^{j,i}=\frac{{\widehat{\omega }}_t^{j,i}}{{\displaystyle {\sum }_{i=1}^N{\widehat{\omega }}_t^{j,i}}}(39)$

Next, the effective sample size, $N_{eff}$ from the weights are calculated, in order to measure the degeneracy problem. If the $N_{eff}$ is smaller than the threshold (e.g.: $N_{eff}$ <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.

$N_{eff}=\frac{1}{{\displaystyle {\sum }_{i=1}^N{\left({\widehat{\omega }}_t^{j,i}\right)}^2}}(40)$

${\omega }_t^{j,i}=1/N(41)$

where ${\widehat{\omega }}_t^{j,i}$ is defined in Eq. 38, ${\omega }_t^{j,i}$ 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

${\widehat{x}}_{t|t}=\frac{1}{M}{\displaystyle \sum _{j=1}^M{\displaystyle \sum _{i=1}^N{\omega }_t^{j,i}{x}_t^{j,i}}}(42)$

where ${\omega }_t^{j,i}$ is defined in Eqs 39, 41, $x_t^{j,i}$ is defined in Eq. 37. For the purpose of next iteration in estimation, the updated measurement estimate, ${\widehat{\zeta }}_{t|t}$ is determined and used to produce the updated measurement error, $e_{\text{W},t|t}^j$ as shown in Eqs 43, 44 respectively.

${\widehat{\zeta }}_{t|t}=g\left({\widehat{x}}_{t|t}\right)(43)$

$e_{\text{W},t|t}^j=W_t^j-{\widehat{\zeta }}_{t|t}(44)$

where $W_t^j$ is defined in Eq. 25, ${\widehat{x}}_{t|t}$ 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 $\left(\Omega ,\mathrm{ℱ},P\right)$ be a probability space where $F=\text{Β}\left({\mathrm{ℝ}}^{n_x}\right)$ is the Borel set of ${\mathrm{ℝ}}^{n_x}$, the Borel set is the standard set of all possible probability events on ${\mathrm{ℝ}}^{n_x}$. Two types of vector-valued stochastic process namely system state, $X=\left\{X_t,t\in \mathrm{\mathbb{N}}\right\}$ and observation, $W=\left\{W_t,t\in \mathrm{\mathbb{N}}\right\}$ are involved in this space. The system state, X is a Markov process of initial distribution X_o∼µ and probability transition kernel, K(x_t|x_t-1).

$p\left(X_t^n\in C_n\text{|}{X}_{t-1}=x_{t-1}\right)={\displaystyle \underset{C_n}{\int }K\left(x_t^n|x_{t-1}\right)d{x}_t^n,C_n\in \left({\mathrm{ℝ}}^{n_x}\right)}\text{for}n=1:n_s(45)$

where $X_t=\left\{X_t^1,X_t^2,\dots ,X_t^{n_s}\right\}$, $X_{t-1}=\left\{X_{t-1}^1,X_{t-1}^2,\dots ,X_{t-1}^{n_s}\right\}$ , n_x is the dimension of the states, n_s is the number of the states. Since there are three states namely flow, stage, and cross section, the current states are represented as $x_t=\left\{Q_t,A_t,H_t\right\}$, and the states at previous time, $x_{t-1}=\left\{Q_{t-1},A_{t-1},H_{t-1}\right\}$. The overall observation is represented as $W_t=W_t^1,\dots ,W_t^{n_k}$ for $1\le k\le n_k$ of n_k sensors that is independent to each other and have marginal distribution

$p\left(W_t^k\in D\text{|}{X}_t=x_t\right)={\displaystyle \underset{D}{\int }p\left(w_t^k|x_t\right)d{w}_t^k,D\in B\left({\mathrm{ℝ}}^{n_w}\right)}(46)$

where $w_t^k$ is the overall observation by kth sensor at time t, $n_w$ is the dimension of the observation and $x_t$ is defined in Eq. 45. For missing observation problem, consider the non-response vector-valued stochastic process, $R=\left\{R_t,t\in \mathrm{\mathbb{N}}\right\}$ with $n_w$ dimensional vector. The availability of the observation is indicated using $r_t^j\in \left(\mathrm{0,1}\right)$ and introduced the following sets

$\begin{array}{l}{\xi }_t=\left\{w_t^k|r_t^k=0for0\le k\le n_w\right\}\\ {\gamma }_t=\left\{w_t^k|r_t^k=1for0\le k\le n_w\right\}\end{array}(47)$

where ${\xi }_t$ is the missing observation at time t, ${\gamma }_t$ is the available observation at time t, $r_t^k$ is the indicator for the availability of the observation, $n_w$ and $w_t^k$ 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

$p\left({\xi }_t\text{|}{r}_{0:t},{\gamma }_{0:t}\right)=p\left({\Xi }_t\in d{\xi }_t|R_{0:t}=r_{0:t},{\Psi }_{0:t}={\gamma }_{0:t}\right)(48)$

where ${\xi }_t$ is the missing observation at time t, ${\gamma }_{0:t}$ is the available observation from time 0 to time t, and $r_{0:t}$ is the indicator for the availability of the observation from time 0 to time t, ${\Xi }_t$ and ${\Psi }_{0:t}$ 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 ${\psi }_t^n$ and the probability density of the states during available observation, ${\eta }_{\beta |\alpha :\beta }^n$ are described as follows:

${\psi }_t^n\triangleq p\left(X_t^n\in d{x}_t^n|{\Psi }_{0:t}={\gamma }_{0:t},R_{0:t}=r_{0:t}\right)(49)$

${\eta }_{\beta |\alpha :\varsigma }^n\triangleq p\left(X_{\beta }^n\in d{x}_{\beta }^n|{\Psi }_{\alpha }={\gamma }_{\alpha },\dots ,{\Psi }_{\varsigma }={\gamma }_{\zeta }\right)(50)$

where ${\Psi }_{\alpha }={\gamma }_{\alpha },\dots ,{\Psi }_{\varsigma }={\gamma }_{\varsigma }$ represent the available observation from time $\alpha$ to time $\varsigma$, $\beta$ is the time for state during the available observation, ${\Psi }_{0:t}$ is defined in Eq. 17, $x_t$ and $X_t$ are defined in Eq. 45, ${\gamma }_{0:t}$, $R_{0:t}$, and $r_{0:t}$ are defined in Eq. 48. For notational convenience, the probability density ${\eta }_{\beta |\alpha :\varsigma }^n$ is written as ${\eta }_{t|t}^n$. So, the distribution of ${\psi }_t^n$ and ${\eta }_{t|t}^n$ are related as follows:

${\psi }_t^n={\displaystyle \int {\eta }_{t|t}^{n}p\left({\xi }_t|{\gamma }_{0:t}\right)d{\xi }_t}(51)$

where ${\eta }_{t|t}^i$ is the probability density of the states at time t, ${\xi }_t,$ and ${\gamma }_{0:t}$ are defined in Eq. 48, ${\psi }_t^n$ is defined in Eq. 49. By applying the standard Bayesian filtering theory, the posterior probability density of the available observation can be expressed as

${\eta }_{t|t}^n=\frac{p\left(w_t|x_t\right){\eta }_{t|t-1}^n}{{\displaystyle \int p\left(w_t|x_t\right){\eta }_{t|t-1}^n}}(52)$

where ${\eta }_{t|t-1}^n$ is the probability density of the states at time t-1, ${\eta }_{t|t}^n$ is defined in Eq. 51, $x_t$ is defined in Eq. 45, $w_t$ 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, ${\eta }_{t|t-1}^n$ and the overall probability density, ${\psi }_t^n$.

${\psi }_t^n={\displaystyle \int \left(\frac{p\left(w_t|x_t\right){\eta }_{t|t-1}^n}{{\displaystyle \int p\left(w_t|x_t\right){\eta }_{t|t-1}^{n}n}}\right)p\left({\xi }_t|{\gamma }_{0:t}\right)d{\xi }_t}(53)$

Since the missing observation is very much affecting the probability distribution, the additional knowledge of this problem is covered by introducing the empirical distribution, ${\left({\eta }_{t|t}^n\right)}^N$ to (51) as a replacement to the true distribution of the probability density, ${\eta }_{t|t}^n$. So, the posterior probability density for each state with N particles is represented as

${\left({\psi }^n\right)}^N={\displaystyle \int {\left({\eta }_{t|t}^n\right)}^{N}p\left({\xi }_t|{\gamma }_t\right)d{\xi }_t}(54)$

where ${\gamma }_t$ is the available observation at time t, ${\left({\eta }_{t|t}^n\right)}^N$ is the probability density of each states during available observation, ${\xi }_t$ is defined in Eq. 48. The empirical distribution, ${\left({\eta }_{t|t}^n\right)}^N$ consists of N particles that are distributed approximately according to ${\eta }_{t|t}^n$.

${\left({\eta }_{t|t}^n\right)}^N={\displaystyle \sum _{i=1}^N{\omega }_t^i\delta {\left(x_t^n\right)}^i}(55)$

given

${\omega }_t^i=\frac{p\left(w_t|{\left(x_t\right)}^i\right)}{{\displaystyle {\sum }_{i=1}^{N}p\left(w_t|{\left(x_t\right)}^i\right)}}(56)$

${\left(x_t\right)}^i$ is the set of particles distributed approximately according to ${\eta }_{t-1|t-1}$ for all state, given by ${\left(x_t\right)}^i\sim K\left(x_t|{\left(x_{t-1}\right)}^i\right)$ as in the standard bootstrap procedure. Next, the missing observation is incorporated into Eq. 32 by applying naive Monte Carlo approximation with ${\xi }_t^j\sim p\left({\xi }_t|{\gamma }_t\right)$ for $0\le j\le M$. The posterior probability density of the desired states with N and M is represented as

${\left({\psi }_t^n\right)}^{N,M}={\left({\eta }_{t|t}^n\right)}^N\frac{1}{M}{\displaystyle \sum _{j=1}^M\delta {\xi }_t^j}(57)$

where ${\left({\eta }_{t|t}^n\right)}^N$ is the posterior probability density for the related states with N particles, ${\xi }_t^j$ 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 $p\left(X_t^n|{\Psi }_{0:t}\right)$ and the true density can be analyzed through convergence analysis.

Part 1: Prediction Using SVSF

Consider the one step prediction error

${\tilde{x}}_{t|t-1}=x_t-{\widehat{x}}_{t|t-1}=f\left(x_{t-1}\right)+{\upsilon }_{t-1}-f\left({\widehat{x}}_{t-1|t-1}\right)(58)$

where ${\widehat{x}}_{t-1|t-1}$ is the estimate of $x_{t-1}$ at time $t-1$ with the initial state, ${\widehat{x}}_{t|t-1}$ is the one step prediction at time $t-1$, $x_t$ and ${\upsilon }_{t-1}$ are defined in Eq. 1.

By using the Taylor series expansion around ${\widehat{x}}_{t-1|t-1}$, $f\left(x_{t-1}\right)$ is linearized as follows

$f\left(x_{t-1}\right)=f\left({\widehat{x}}_{t-1|t-1}\right)+A_{t-1}{\tilde{x}}_{t-1|t-1}+\varphi {\tilde{x}}_{t-1|t-1}^2(59)$

where $A_{t-1}=\frac{\partial f\left(x_{t-1}\right)}{\partial x_{t-1}}{|}_{x_{t-1}={\widehat{x}}_{t-1|t-1}}$, and $\varphi {\tilde{x}}_{t-1|t-1}^2$ is the high order terms of the Taylor series expansion

$\varphi {\tilde{x}}_{t-1|t-1}^2={\text{Λ}}_{t-1}{\aleph }_{1,t-1}{L}_{t-1}{\tilde{x}}_{t-1|t-1}(60)$

where ${\Lambda }_{t-1}$ is a problem dependent scaling matrix, $L_{t-1}$ is introduce to provide extra degree of freedom to the filter, ${\aleph }_{1,t-1}$ is an unknown time-varying matrix for linearization errors of the dynamical model that satisfies

${\aleph }_{1,t-1}{\aleph }_{1,t-1}^T\le I(61)$

that gives

${\tilde{x}}_{t|t-1}=\left(A_{t-1}+{\text{Λ}}_{t-1}{\aleph }_{1,t-1}{L}_{t-1}\right){\tilde{x}}_{t-1|t-1}+{\upsilon }_{t-1}(62)$

where $A_{t-1}$ is defined in Eq. 59, ${\upsilon }_{t-1}$ is defined in Eq. 1, ${\Lambda }_{t-1}$ is defined in Eq. 60, ${\aleph }_{1,t-1}$ and $L_{t-1}$ are defined in Eq. 60. Applying the Taylor series expansion around ${\widehat{x}}_{t|t-1}$ to the $g\left(x_{t-1}\right)$, the innovation of the filter is represented as

${\tilde{W}}_t={\text{W}}_t-g\left({\widehat{x}}_{t|t-1}\right)=\left(B_t+{\text{Ω}}_{\text{t}}{\aleph }_{2,t}{L}_t\right){\tilde{x}}_{t|t-1}+{\epsilon }_t(63)$

where $B_t=\frac{\partial g\left(x_t\right)}{\partial x_t}{|}_{x_t={\widehat{x}}_{t|t}}$, ${\Omega }_t$ is a problem dependent scaling matrix, $L_t$ is introduce to provide extra degree of freedom to the filter, ${\aleph }_{2,t}$ is an unknown time-varying matrix for linearization errors of the dynamical model that satisfies ${\aleph }_{2,t}{\aleph }_{2,t}^T\le I$ as in Eq. 61, ${\tilde{x}}_{t|t-1}$ is defined in Eq. 62, ${\widehat{x}}_{t|t-1}$ is defined in Eq. 59, ${\text{W}}_t$ and ${\epsilon }_t$ are defined in Eq. 2. The filtering error is represented as

${\tilde{x}}_{t|t}=x_t-{\widehat{x}}_{t|t}=\left(I-K_{svsf,t}\left(B_t+{\text{Ω}}_{\text{t}}{\aleph }_{2,t}{L}_t\right)\right){\tilde{x}}_{t|t-1}-K_{svsf,t}{\epsilon }_t(64)$

where $K_{svsf,t}$ is the gain, ${\widehat{x}}_{t|t}$ is the estimated state from SVSF, $x_t$ is defined in Eq. 45, ${\widehat{x}}_{t|t-1}$ is defined in Eq. 59, $B_t, {\Omega }_{t, }{\aleph }_{2,t}, L_t$ are defined in Eq. 63, ${\epsilon }_t$ is defined in Eq. 2,

Theorem 1. The one step prediction error covariance, $Cov{P}_{t|t-1}$ is given by

$Cov{P}_{t|t-1}=\left(A_{t-1}+{\text{Λ}}_{t-1}{\aleph }_{1,t-1}{L}_{t-1}\right)Cov{P}_{t-1|t-1}{\left(A_{t-1}+{\text{Λ}}_{t-1}{\aleph }_{1,t-1}{L}_{t-1}\right)}^T+Cov{Q}_{t-1}(65)$

where $Cov{P}_{t-1|t-1}$ is the previous prediction covariance, $Cov{Q}_{t-1}$ is defined in Eq. 23, $A_{t-1}$ is defined in Eq. 59, ${\upsilon }_{t-1}$ is defined in Eq. 1, ${\Lambda }_{t-1}$ is defined in Eq. 60, ${\aleph }_{1,t-1},$ and $L_{t-1}$ are defined in Eq. 60

Theorem 2. The filtering error covariance $Cov{P}_{t|t}$ is given by

$Cov{P}_{t|t}=\left(I-\left(K_{svsf,t}\left(B_t+{\text{Ω}}_{\text{t}}{\aleph }_{2,t}{L}_t\right)\right)\right)Cov{P}_{t|t-1}\times {\left(I-\left(K_{svsf,t}\left(B_t+{\text{Ω}}_{\text{t}}{\aleph }_{2,t}{L}_t\right)\right)\right)}^T+\left(K_{svsf,t}\left({\text{CovR}}_{\text{t}}\right)K_{svsf,t}^T\right)(66)$

where $\mathbb{E}\left\{{\epsilon }_t{\epsilon }_t^T\right\}=Cov{R}_t$ is defined in Eq. 23, $K_{svsf,t}$ is defined in Eq. 64, $B_t, {\Omega }_{t, }{\aleph }_{2,t}, L_t$ are defined in Eq. 63.

Lemma 1. Given matrices A, H, E, and F with appropriate dimensions such that $F{F}^T\le I$. Let X be a symmetric positive definite matrix and $\gamma$ be an arbitrary positive constant such that ${\gamma }^{-1}I-EX{E}^T>0$. Then the following inequality holds

$\left(A+HFE\right)X{\left(A+HFE\right)}^T\le A\left(X^{-1}-\gamma E^{T}E\right)A^T+{\gamma }^{-1}H{H}^T(67)$

Theorem 3: Consider Theorem 1 and Theorem 2 and assume ${\aleph }_{1,t-1}{\aleph }_{1,t-1}^T\le I$ and ${\aleph }_{2,t}{\aleph }_{2,t}^T\le I$ as in Eqs 60, 62 are true. Let ${\gamma }_{1,t-1},{\gamma }_{2,t}$ 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:

${\text{Σ}}_{t|t-1}=A_{t-1}{\left({\text{Σ}}_{t-1|t-1}^{-1}-{\gamma }_{1,t-1}^{-1}{L}_{t-1}{L}_{t-1}^T\right)}^{-1}{A}_{t-1}^T+{\gamma }_{1,t-1}^{-1}{\text{Λ}}_{t-1}{\text{Λ}}_{t-1}^T+Cov{Q}_{t-1}(68)$

${\text{Σ}}_{t|t}=\left(I-\left(K_{s,t}{B}_t\right)\right){\left({\displaystyle {\sum }_{t|t-1}^{-1}-{\gamma }_{2,t}^{-1}{L}_t{L}_t^T}\right)}^{-1}{\left(I-\left(K_{s,t}{B}_t\right)\right)}^T+{\gamma }_{2,t}^{-1}{K}_{s,t}{\Omega }_{\text{t}}{\Omega }_{\text{t}}^T{K}_{s,t}^T+\left(K_{s,t}\left({\text{CovR}}_{\text{t}}\right)K_{s,t}^T\right)(69)$

where $K_{s,t}$ is the gain like $K_{svsf,t}$ in Eq. 64, $Cov{Q}_{t-1}$ and $Cov{R}_t$ are defined in Eq. 23, $A_{t-1}$ is defined in Eq. 59, ${\Lambda }_{t-1}$ is defined in Eq. 60, $L_{t-1}$ is defined in Eq. 60, $B_t, {\Omega }_{\text{t},}{\aleph }_{2,t}, L_t$ are defined in Eq. 63, ${\gamma }_{1,t-1}, {\gamma }_{2,t}$ are positive scalars.With initial condition ${\text{Σ}}_{0|0}=Cov{P}_{0|0}>0$ have positive definite solutions ${\text{Σ}}_{t|t-1}$ and ${\text{Σ}}_{t|t}$ such that for all $0<t<N_t$ the following two constraints

${\gamma }_{1,t-1}^{-1}I-L_{t-1}{\text{Σ}}_{t-1|t-1}{L}_{t-1}^T>0(70)$

${\gamma }_{2,t}^{-1}I-L_t{\text{Σ}}_{t|t-1}{L}_t^T>0(71)$

are satisfied.The filter gain for the upper bound $K_{s,t}$ is referring to the $K_{svsf,t}$ in Eq. 32 and minimizes the upper bound so that

$Cov{P}_{t|t}\le {\text{Σ}}_{t|t}(72)$

Lemma 2. For $0<k<N$, suppose that $X=X^T>0,S_k\left(X\right)=S_k^T\left(X\right)\in {\mathrm{ℝ}}^{nxn}$ and $H_k\left(X\right)=H_k^T\left(X\right)\in {\mathrm{ℝ}}^{nxn}$. If

$S_k\left(Y\right)\ge S_k\left(X\right),\forall X\le Y=Y^T(73)$

and

$H_k\left(Y\right)\ge S_k\left(Y\right)(74)$

Then the solutions for $M_{k}and{N}_k$ to the following difference equations are

$M_{k+1}=S_k\left(M_k\right), N_{k+1}=H_k\left(N_k\right), M_o=N_0>0(75)$

Satisfy. $M_k\le N_k$

Proof. By referring to Lemma 2, rewrite the error covariance matrices as the function of $Cov{P}_{t|t-1}$.

$Cov{P}_{t|t-1}\left(Cov{P}_{t-1|t-1}\right)=\left(A_{t-1}+{\text{Λ}}_{t-1}{\aleph }_{1,t-1}{L}_{t-1}\right)Cov{P}_{t-1|t-1}{\left(A_{t-1}+{\text{Λ}}_{t-1}{\aleph }_{1,t-1}{L}_{t-1}\right)}^T+Cov{Q}_{t-1}(76)$

$Cov{P}_{t|t}\left(Cov{P}_{t|t-1}\right)=\left(I-\left(K_{svsf,t}\left(B_t+{\text{Ω}}_{\text{t}}{\aleph }_{2,t}{L}_t\right)\right)\right).Cov{P}_{t|t-1} \times {\left(I-\left(K_{svsf,t}\left(B_t+{\text{Ω}}_{\text{t}}{\aleph }_{2,t}{L}_t\right)\right)\right)}^T+\left(K_{svsf,t}\left({\text{CovR}}_{\text{t}}\right)K_{svsf,t}^T\right)(77)$

where $Cov{P}_{t-1|t-1}$ is defined in Eq. 65, $Cov{Q}_{t-1}$ , $Cov{R}_t$ are defined in Eq. 23, $A_{t-1}$ is defined in Eq. 59, ${\upsilon }_{t-1}$ is defined in Eq. 1), ${\Lambda }_{t-1}$, ${\aleph }_{1,t-1}$, and $L_{t-1}$ are defined in Eq. 60, $K_{svsf,t}$ is defined in Eq. 56, $B_t, {\Omega }_{\text{t},}{\aleph }_{2,t}, L_t$ are defined in Eq. 63.Next, rewrite the upper bound obtained in Eq. 67 and Eq. 68 as the function of ${\text{Σ}}_{t-1|t-1}$ and ${\text{Σ}}_{t|t-1}$

${\text{Σ}}_{t|t-1}\left({\text{Σ}}_{t-1|t-1}\right)=A_{t-1}{\left({\text{Σ}}_{t-1|t-1}^{-1}-{\gamma }_{1,t-1}^{-1}{L}_{t-1}{L}_{t-1}^T\right)}^{-1}{A}_{t-1}^T+{\gamma }_{1,t-1}^{-1}{\text{Λ}}_{t-1}{\text{Λ}}_{t-1}^T+Cov{Q}_{t-1}(78)$

${\text{Σ}}_{t|t}\left({\text{Σ}}_{t|t-1}\right)=\left(I-\left(K_{s,t}{B}_t\right)\right){\left({\text{Σ}}_{t|t-1}^{-1}-{\gamma }_{2,t}^{-1}{L}_t{L}_t^T\right)}^{-1}{\left(I-\left(K_{s,t}{B}_t\right)\right)}^T+{\gamma }_{2,t}^{-1}{K}_{s,t}{\Omega }_{\text{t}}{\Omega }_{\text{t}}^T{K}_{s,t}^T+\left(K_{s,t}\left({\text{CovR}}_{\text{t}}\right)K_{s,t}^T\right)(79)$

where $K_{s,t}$ is defined in Eq. 61, $Cov{Q}_{t-1}$ and $Cov{R}_t$ are defined in Eq. 23, $A_{t-1}$ is defined in Eq. 59, ${\Lambda }_{t-1}$, $L_{t-1}$ is defined in Eq. 60, $B_t, {\Omega }_{\text{t},}{\aleph }_{2,t}, L_t$ are defined in Eq. 63, ${\gamma }_{1,t-1},{\gamma }_{2,t}$ are positive scalars.Consider Eqs 75–78, the error covariance matrices and the upper bound satisfy condition Eq. 72, 73 in Lemma 2.

$Cov{P}_{t|t}\left(Cov{P}_{t|t-1}\right)\ge Cov{P}_{t|t-1}\left(Cov{P}_{t-1|t-1}\right)(80)$

${\displaystyle {\sum }_{t|t}\left({\text{Σ}}_{t|t-1}\right)}\ge Cov{P}_{t|t}\left(Cov{P}_{t|t-1}\right)(81)$

which gives $Cov{P}_{t|t}\le {\text{Σ}}_{t|t}$.

Part 2: Prediction Using MIPF

Consider: if ${\left({\mu }_N\right)}_{N=1}^{\infty }$ is a sequence of random probability measures, then ${\mu }_N$ converges to $\mu \in B\left({\mathrm{ℝ}}^{n_x}\right)$ if for any continuous bounded function $\phi \in B\left({\mathrm{ℝ}}^{n_x}\right)$

$\underset{N\to \infty }{\lim }\mathbb{E}\left[{\left(\left({\mu }_N,\phi \right)-\left(\mu ,\phi \right)\right)}^2\right]=0$

Assumption. The likelihood function $p\left(W_t\text{|}{X}_t\right)$ is a bounded function in argument and represented as $\pi <\infty$.

Lemma 3. For any $\phi \in B\left({\mathrm{ℝ}}^{n_x}\right)$ with random variables ${\left\{x_t^i\right\}}_{i=1}^N$ obtained from states, ${\widehat{x}}_{t|t}$ and prediction error covariance matrices, $Cov{P}_{t|t}$ from SVSF as in Eq. 37.

$\begin{array}{l}\mathbb{E}{\left|\left({\left({\eta }_{t-1|t-1}^n\right)}^N,\phi \right)-\left({\eta }_{t-1|t-1}^n,\phi \right)\right|}^4=\frac{1}{N^4}\mathbb{E}{\left[{\displaystyle \sum _{i=1}^N\left(f\left(x_t^i\right)-\mathbb{E}\left(f\left(x_t^i\right)\right)\right)}\right]}^4 \le \frac{2}{N^2}\mathbb{E}{\left[f\left(x_t^i\right)-\mathbb{E}\left(f\left(x_t^i\right)\right)\right]}^4\\ \le C_{t-1|t-1}\frac{{\Vert \phi \Vert }^2}{N}forn=1:n_s\end{array}(82)$

Then, for any $\phi \in B\left({\mathrm{ℝ}}^{n_x}\right)$

$\mathbb{E}{\left|\left({\left({\eta }_{t|t-1}^n\right)}^N,\phi \right)-\left({\eta }_{t|t-1}^n,\phi \right)\right|}^2=C_{t|t-1}^n\frac{{\Vert \phi \Vert }^2}{N}(83)$

where ${\eta }_{t-1|t-1}^n$ is the probability density of the random variables, ${\eta }_{t|t-1}^n$ is the probability density of the one step state estimate, $C_{t|t-1}^n$ and $C_{t-1|t-1}^n$ are the constant, $x_t^i$ is the particles obtained, and$n_s$ is the number of the states; N is described in Eq. 37.

ProofConsider

$\begin{array}{l}\left|\left({\left({\eta }_{t|t-1}^n\right)}^N,\phi \right)-\left({\eta }_{t|t-1}^n,\phi \right)\right|\\ \le \left|\left({\left({\eta }_{t|t-1}^n\right)}^N,\phi \right)-\left({\left({\eta }_{t-1|t-1}^n\right)}^N,K\phi \right)\right|+\left|\left({\left({\eta }_{t-1|t-1}^n\right)}^N,K\phi \right)-\left({\eta }_{t-1|t-1}^n,K\phi \right)\right|\end{array}(84)$

Let ${\mathcal{G}}_{t-1}$ be the $\sigma -field$ generated by ${\left\{x_{t-1}^i\right\}}_{i=1}^N$, then

$\mathbb{E}\left[\left({\eta }_{t|t-1}^N,\phi \right)|{\mathcal{G}}_{t-1}\right]=\left({\eta }_{t-1|t-1}^N,K\phi \right)(85)$

and, as $K\phi \le \phi$

$\mathbb{E}\left[{\left(\left({\left({\eta }_{t|t-1}^n\right)}^N,\phi \right)-E\left(\left({\left({\eta }_{t|t-1}^n\right)}^N,K\phi \right)\right)|G_{t-1}\right)}^2|G_{t-1}\right]\le \frac{{\Vert \phi \Vert }^2}{N}(86)$

Using Minkowski’s inequality

$\mathbb{E}{\left[{\left(\left({\left({\eta }_{t|t-1}^n\right)}^N,\phi \right)-\left({\eta }_{t|t-1}^n,\phi \right)\right)}^2\right]}^{\frac{1}{2}}\le \sqrt{C_{t-1|t-1}^n}\frac{\phi }{\sqrt{N}}+\frac{\phi }{\sqrt{N}}\le C_{t|t-1}^n\frac{\Vert \phi \Vert }{\sqrt{N}}(87)$

where $C_{t|t-1}^n={\left(\sqrt{C_{t-1|t-1}^n}+1\right)}^2$.

Part 3: Updating Using MIPF

Lemma 4. For any $\phi \in B\left({ℝ}^{n_x}\right)$

$\mathbb{E}{\left[\left({\left({\eta }_{t|t-1}^n\right)}^N,\phi \right)-\left({\eta }_{t|t-1}^n,\phi \right)\right]}^2\le C_{t|t-1}^n\frac{{\Vert \phi \Vert }^2}{N} (88)$

$where{C}_{t|t-1}^n={\left(\sqrt{C_{t-1|t-1}^n}+1\right)}^2$Then for any $\phi \in B\left({\mathrm{ℝ}}^{n_x}\right)$

$\mathbb{E}{\left[\left({\left({\tilde{\eta }}_{t|t}^n\right)}^N,\phi \right)-\left({\eta }_{t|t}^n,\phi \right)\right]}^2={\tilde{C}}_{t|t-1}^n\frac{{\Vert \phi \Vert }^2}{N}(89)$

where ${\eta }_{t|t-1}^n$ is the probability density of the one step state estimate, ${\tilde{\eta }}_{t|t}^n$ is the probability density of the state without consideration of missing data, $C_{t|t-1}^n$ and ${\tilde{C}}_{t|t-1}^n$ are the constant, N is described in Eq. 37.

Proof: Consider

$\left({\left({\tilde{\eta }}_{t|t}^n\right)}^N,\phi \right)-\left({\eta }_{t|t}^n,\phi \right)=\frac{\left({\left({\eta }_{t|t-1}^n\right)}^N,\pi \phi \right)}{\left({\left({\eta }_{t|t-1}^n\right)}^N,\pi \right)}-\frac{\left({\eta }_{t|t-1}^n,\pi \phi \right)}{\left({\eta }_{t|t-1}^n,\pi \right)}(90)$

For the first part of Eq. 90

$\left|\frac{\left({\left({\eta }_{t|t-1}^n\right)}^N,\pi \phi \right)}{\left({\left({\eta }_{t|t-1}^n\right)}^N,\pi \right)}-\frac{\left({\left({\eta }_{t|t-1}^n\right)}^N,\pi \phi \right)}{\left({\eta }_{t|t-1}^n,\pi \right)}\right|\le \frac{\phi }{\left({\eta }_{t|t-1}^n,\pi \right)}\left|\left({\eta }_{t|t-1}^n,\pi \right)-\left({\left({\eta }_{t|t-1}^n\right)}^N,\pi \right)\right|(91)$

Using Minkowski’s inequality

$\mathbb{E}{\left[{\left(\left({\left({\tilde{\eta }}_{t|t}^n\right)}^N,\phi \right)-\left({\eta }_{t|t}^n,\phi \right)\right)}^2\right]}^{\frac{1}{2}}\le \frac{\Vert \phi \Vert }{\left({\left({\eta }_{t|t-1}^n\right)}^N,\pi \right)}\mathbb{E}{\left[{\left(\left({\eta }_{t|t-1}^n,\pi \phi \right)-\left({\left({\eta }_{t|t-1}^n\right)}^N,\pi \right)\right)}^2\right]}^{\frac{1}{2}}+\frac{\mathbb{E}{\left[{\left(\left({\left({\eta }_{t|t-1}^n\right)}^N,\pi \phi \right)-\left({\eta }_{t|t-1}^n,\pi \phi \right)\right)}^2\right]}^{\frac{1}{2}}}{\left({\left({\eta }_{t|t-1}^n\right)}^N,\pi \right)}\le \frac{2\sqrt{C_{t|t-1}^n}\Vert \pi \Vert }{\left({\left({\eta }_{t|t-1}^n\right)}^N,\pi \right)}\frac{\Vert \phi \Vert }{\sqrt{N}}(92)$

Part 4: Updating With New Imputation Using MIPF

Lemma 5. For any $\phi \in B\left({\mathrm{ℝ}}^{n_x}\right)$

$\mathbb{E}{\left[\left({\left({\tilde{\eta }}_{t|t}^n\right)}^N,\phi \right)-\left({\eta }_{t|t}^n,\phi \right)\right]}^2\le {\tilde{C}}_{t|t-1}^n\frac{{\Vert \phi \Vert }^2}{N} (93)$

$\text{where} {\tilde{C}}_{t|t-1}^n= {\left(\frac{2\sqrt{C_{t|t-1}^n}\pi }{\left({\left({\eta }_{t|t-1}^n\right)}^N,\pi \right)}\right)}^2$. Then, for any $\phi \in B\left({\mathrm{ℝ}}^{n_x}\right)$ that includes new imputation $\mathit{\varrho }\in B\left({\mathrm{ℝ}}^{n_z}\right)$

$\mathbb{E}{\left[\left({\left({\widetilde{\text{Ψ}}}_{t|t}^n\right)}^{N,M},\phi \right)-\left({\text{Ψ}}_{t|t}^n,\phi \right)\right]}^2={\tilde{f}}_{t|t}^n\frac{{\Vert \phi \Vert }^2}{MN}(94)$

where ${\tilde{\eta }}_{t|t}^n$ is the probability density of the state without the consideration of missing data, ${\tilde{\Psi }}_{t|t}^n$ is the probability density of the state with the consideration of missing data, and${\tilde{C}}_{t|t-1}^n$ and ${\tilde{f}}_{t|t}^n$ are the constant; M and N are described in Eq. 37.

Proof: Consider

$\left({\left({\tilde{\eta }}_{t|t}^n\right)}^N{H}_t^M,\phi \right)-\left({\eta }_{t|t}^n{H}_t,\phi \right)=\frac{\left({\left({\eta }_{t|t-1}^n\right)}^N{H}_t^M,\mathit{\varrho }\pi \phi \right)}{\left({\left({\eta }_{t|t-1}^n\right)}^N{H}_t^M,\mathit{\varrho }\pi \right)}-\frac{\left({\left({\eta }_{t|t-1}^n\right)}^N{H}_t^M,\mathit{\varrho }\pi \phi \right)}{\left({\eta }_{t|t-1}^n{H}_t,\mathit{\varrho }\pi \right)}+\frac{\left({\left({\eta }_{t|t-1}^n\right)}^N{H}_t^M,\mathit{\varrho }\pi \phi \right)}{\left({\eta }_{t|t-1}^n{H}_t,\mathit{\varrho }\pi \right)}-\frac{\left({\eta }_{t|t-1}^n{H}_t,\mathit{\varrho }\pi \phi \right)}{\left({\eta }_{t|t-1}^n{H}_t,\mathit{\varrho }\pi \right)} (95)$

For the first part of Eq. 95

$\left|\frac{\left({\left({\eta }_{t|t-1}^n\right)}^N{H}_t^M,\mathit{\varrho }\pi \phi \right)}{\left({\left({\eta }_{t|t-1}^n\right)}^N{H}_t^M,\mathit{\varrho }\pi \right)}-\frac{\left({\left({\eta }_{t|t-1}^n\right)}^N{H}_t^M,\mathit{\varrho }\pi \phi \right)}{\left({\eta }_{t|t-1}^n{H}_t,\mathit{\varrho }\pi \right)}\right|\le \frac{\Vert \phi \Vert }{\left({\eta }_{t|t-1}^n{H}_t,\mathit{\varrho }\pi \right)}\left|\left({\eta }_{t|t-1}^n{H}_t,\mathit{\varrho }\pi \right)-\left({\left({\eta }_{t|t-1}^n\right)}^N{H}_t^M,\mathit{\varrho }\pi \right)\right|(96)$

For the second part of Eq. 94 that refer to Eq. 89

${\left|\frac{\left({\left({\eta }_{t|t-1}^n\right)}^N{H}_t^M,\mathit{\varrho }\pi \phi \right)}{\left({\eta }_{t|t-1}^n{H}_t,\mathit{\varrho }\pi \right)}-\frac{\left({\eta }_{t|t-1}^n{H}_t,\mathit{\varrho }\pi \phi \right)}{\left({\eta }_{t|t-1}^n{H}_t,\mathit{\varrho }\pi \right)}\right|}^2={\tilde{f}}_{t|t-1}^n\frac{{\Vert \phi \Vert }^2}{MN}(97)$

Using Minkowski’s inequality

$\mathbb{E}{\left[{\left(\left({\left({\tilde{\eta }}_{t|t}^n\right)}^N{H}_t^M,\phi \right)-\left({\eta }_{t|t}^n{H}_t,\phi \right)\right)}^2\right]}^{\frac{1}{2}}\le \frac{2\sqrt{{\tilde{f}}_{t|t-1}^n}\mathit{\varrho }\pi }{\left({\eta }_{t|t-1}^n{H}_t,\mathit{\varrho }\pi \right)} \frac{\Vert \phi \Vert }{\sqrt{MN}}(98)$

Part 5: Resampling Using MIPF

Lemma 6. For any $\phi \in B\left({\mathrm{ℝ}}^{n_x}\right)$ and $\mathit{\varrho }\in B\left({\mathrm{ℝ}}^{n_z}\right)$

$\mathbb{E}{\left|\left({\left({\widetilde{\text{Ψ}}}_{t|t}^n\right)}^{N,M},\phi \right)-\left({\text{Ψ}}_{t|t}^n,\phi \right)\right|}^2\le {\tilde{f}}_{t|t}^n\frac{{\Vert \phi \Vert }^2}{MN}with{\tilde{f}}_{t|t}^n={\left(\frac{2\sqrt{{\tilde{f}}_{t|t-1}^n}\mathit{\varrho }\pi }{\left({\eta }_{t|t-1}^n{H}_t,\mathit{\varrho }\pi \right)}\right)}^2(99)$

There exists a constant $f_{t|t}^n$, such that for any $\phi \in B\left({\mathrm{ℝ}}^{n_x}\right)$ and $\mathit{\varrho }\in B\left({\mathrm{ℝ}}^{n_z}\right)$

$\mathbb{E}{\left|\left({\left({\text{Ψ}}_{t|t}^n\right)}^{N,M},\phi \right)-\left({\text{Ψ}}_{t|t}^n,\phi \right)\right|}^2=f_{t|t}^n\frac{{\Vert \phi \Vert }^2}{MN}(100)$

where ${\tilde{\Psi }}_{t|t}^n$ is the probability density of the state with the consideration of missing data, ${\Psi }_{t|t}^n$ is the posterior probability density of the states, and ${\tilde{C}}_{t|t-1}^n$ and ${\tilde{f}}_{t|t}^n$ are the constant; M and N are described in Eq. 37.

ProofConsider the following

$\left({\left({\text{Ψ}}_{t|t}^n\right)}^{N,M},\phi \right)-\left({\text{Ψ}}_{t|t}^n,\phi \right)=\left({\left({\text{Ψ}}_{t|t}^n\right)}^{N,M},\phi \right)-\left({\left({\widetilde{\text{Ψ}}}_{t|t}^n\right)}^{N,M},\phi \right)+\left({\left({\widetilde{\text{Ψ}}}_{t|t}^n\right)}^{N,M},\phi \right)-\left({\text{Ψ}}_{t|t}^n,\phi \right)(101)$

By using Minkowski’s inequality

$\begin{array}{l}\mathbb{E}{\left[{\left(\left({\left({\text{Ψ}}_{t|t}^n\right)}^{N,M},\phi \right)-\left({\text{Ψ}}_{t|t}^n,\phi \right)\right)}^2\right]}^{\frac{1}{2}}\\ \le \mathbb{E}{\left[{\left(\left({\left({\text{Ψ}}_{t|t}^n\right)}^{N,M},\phi \right)-\left({\left({\widetilde{\text{Ψ}}}_{t|t}^n\right)}^{N,M},\phi \right)\right)}^2\right]}^{\frac{1}{2}}\mathrm{ }+\mathbb{E}{\left[{\left(\left({\left({\widetilde{\text{Ψ}}}_{t|t}^n\right)}^{N,M},\phi \right)-\left({\text{Ψ}}_{t|t}^n,\phi \right)\right)}^2\right]}^{\frac{1}{2}}\end{array}(102)$

Let $F_{t-1}$ be the $\sigma -field$ generated by ${\left\{{\tilde{x}}_{t-1}^i\right\}}_{i=1}^N$, then

$\mathbb{E}\left[\left({\left({\text{Ψ}}_{t|t}^n\right)}^{N,M},\phi \right)|F_{t-1}\right]=\left({\left({\widetilde{\text{Ψ}}}_{t|t}^n\right)}^{N,M},\phi \right)(103)$

and

$\mathbb{E}\left[{\left(\left({\left({\text{Ψ}}_{t|t}^n\right)}^{N,M},\phi \right)-\left({\left({\widetilde{\text{Ψ}}}_{t|t}^n\right)}^{N,M},\phi \right)\right)}^2|F_{t-1}\right]=\frac{F}{MN}{\Vert \phi \Vert }^2(104)$

That gives

$\mathbb{E}{\left[{\left(\left({\left({\text{Ψ}}_{t|t}^n\right)}^{N,M},\phi \right)-\left({\text{Ψ}}_{t|t}^n,\phi \right)\right)}^2\right]}^{\frac{1}{2}}\le \frac{\sqrt{F}+{\tilde{f}}_{t|t}^n}{\sqrt{MN}}\Vert \phi \Vert (105)$

Theorem 4. For all $t>0$ there exist $f_{t|t}^n$ that is independent of N but being influenced by M for any bounded function, $\phi \in B\left({\mathrm{ℝ}}^{n_x}\right)$

$\mathbb{E}\left[{\left(\left({\left({\text{Ψ}}_{t|t}^n\right)}^{N,M},\phi \right)-\left({\text{Ψ}}_{t|t}^n,\phi \right)\right)}^2\right]\le f_{t|t}^n\frac{{\Vert \phi \Vert }^2}{MN}with f_{t|t}^n={\left(\sqrt{F}+{\tilde{f}}_{t|t}^n \right)}^2(106)$

$\mathbb{E}\left[{\left(\left({\left({\text{Ψ}}_{t|t}^n\right)}^{N,M},\phi \right)-\left({\text{Ψ}}_{t|t}^n,\phi \right)\right)}^2\right]\le \frac{1}{M^2{N}^2}\mathbb{E}{\left[{\displaystyle \sum _{j=1}^M{\displaystyle \sum _{i=1}^N\left(f\left(x_t^{j,i}\right)-\mathbb{E}\left(f\left(x_t^{j,i}\right)\right)\right)}}\right]}^2(107)$

where ${\Psi }_{t|t}^n$ is the posterior probability density of the states, $f_{t|t}^n$ is the constant, M and N are described in Eq. 37, $x_t^{j,i}$ 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, $n_z$ and dimension of the states, $n_x$ 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

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

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

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

TABLE 2

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

TABLE 3

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

TABLE 4

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

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.

Data Availability Statement

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

Author Contributions

ZI and NJ: Conceptualization; Methodology; Validation; Investigation; Visualization; Writing—original draft. ZI: Formal analysis; Methodology; Validation; Investigation; Visualization; Writing—original draft.

Funding

The research has been carried out under program Research Excellence Consortium (JPT (BPKI) 1000/016/018/25 (57)) provided by Ministry of Higher Education Malaysia (MOHE).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

The author also would like to thank Universiti Teknologi Malaysia (UTM) for work and facilities support.

References

Abaza, M., Anctil, F., Fortin, V., and Turcotte, R. (2014). Sequential Streamflow Assimilation for Short-Term Hydrological Ensemble Forecasting. J. Hydrol. 519, 2692–2706. doi:10.1016/j.jhydrol.2014.08.038

CrossRef Full Text | Google Scholar

Adnan, R. M., Yuan, X., Kisi, O., and Anam, R. (2017). Improving Accuracy of River Flow Forecasting Using LSSVR with Gravitational Search Algorithm. Adv. Meteorol. 2017. doi:10.1155/2017/2391621

CrossRef Full Text | Google Scholar

Arulampalam, M. S., Maskell, S., Gordon, N., and Clapp, T. (2002). A Tutorial on Particle Filters for Online Nonlinear/non-Gaussian Bayesian Tracking. IEEE Trans. Signal. Process. 50 (2), 174–188. doi:10.1109/78.978374

CrossRef Full Text | Google Scholar

Aswathy, S., Sajikumar, N., and Mehsa, M. (2016). Watershed Modelling Using Control System Concept. Proced. Tech. 24, 39–46. doi:10.1016/j.protcy.2016.05.007

CrossRef Full Text | Google Scholar

Bertino, L., Evensen, G., and Wackernagel, H. (2003). Sequential Data Assimilation Techniques in Oceanography. Int. Stat. Rev. 71 (2), 223–241. doi:10.1111/j.1751-5823.2003.tb00194.x

CrossRef Full Text | Google Scholar

Bjerklie, D. M., Birkett, C. M., Jones, J. W., Carabajal, C., Rover, J. A., Fulton, J. W., et al. (2018). Satellite Remote Sensing Estimation of River Discharge: Application to the Yukon River Alaska. J. Hydrol. 561 (April), 1000–1018. doi:10.1016/j.jhydrol.2018.04.005

CrossRef Full Text | Google Scholar

Blythe, T. L., and Schmidt, J. C. (2018). Estimating the Natural Flow Regime of Rivers with Long-Standing Development: The Northern Branch of the Rio Grande. Water Resour. Res. 54 (2), 1212–1236. doi:10.1002/2017wr021919

CrossRef Full Text | Google Scholar

Cao, Y., Ye, Y., Liang, L., Zhao, H., Jiang, Y., Wang, H., et al. (2019). A Modified Particle Filter‐Based Data Assimilation Method for a High‐Precision 2‐D Hydrodynamic Model Considering Spatial‐temporal Variability of Roughness: Simulation of Dam‐Break Flood Inundation. Water Resour. Res. 55, 6049–6068. doi:10.1029/2018wr023568

CrossRef Full Text | Google Scholar

Chai, T., and Draxler, R. R. (2014). Root Mean Square Error (RMSE) or Mean Absolute Error (MAE)? - Arguments against Avoiding RMSE in the Literature. Geosci. Model. Dev. 7 (3), 1247–1250. doi:10.5194/gmd-7-1247-2014

CrossRef Full Text | Google Scholar

Crisan, D., and Doucet, A. (2002). A Survey of Convergence Results on Particle Filtering Methods for Practitioners. IEEE Trans. Signal. Process. 50 (3), 736–746. doi:10.1109/78.984773

CrossRef Full Text | Google Scholar

Ding, D., Wang, Z., Alsaadi, F. E., and Shen, B. (2015). Receding Horizon Filtering for a Class of Discrete Time-Varying Nonlinear Systems with Multiple Missing Measurements. Int. J. Gen. Syst. 44 (2), 198–211. doi:10.1080/03081079.2014.973732

CrossRef Full Text | Google Scholar

Feng, J., Wang, Z., and Zeng, M. (2011). Recursive Robust Filtering with Finite-step Correlated Process Noises and Missing Measurements. Circuits Syst. Signal. Process. 30 (6), 1355–1368. doi:10.1007/s00034-011-9289-6

CrossRef Full Text | Google Scholar

Gadsden, S. A., Al-Shabi, M., Arasaratnam, I., and Habibi, S. R. (2014). Combined Cubature Kalman and Smooth Variable Structure Filtering: A Robust Nonlinear Estimation Strategy. Signal. Process. 96 (PART B), 290–299. doi:10.1016/j.sigpro.2013.08.015

CrossRef Full Text | Google Scholar

Gadsden, S. A., Habibi, S. R., and Kirubarajan, T. (2012). The Smooth Particle Variable Structure Filter. Trans. Can. Soc. Mech. Eng. 36 (2), 177–193. doi:10.1139/tcsme-2012-0013

CrossRef Full Text | Google Scholar

Habibi, S. (2007). The Smooth Variable Structure Filter. Proc. IEEE 95 (5), 1026–1059. doi:10.1109/jproc.2007.893255

CrossRef Full Text | Google Scholar

He, X., Sithiravel, R., Tharmarasa, R., Balaji, B., and Kirubarajan, T. (2014). A Spline Filter for Multidimensional Nonlinear State Estimation. Signal. Process. 102, 282–295. doi:10.1016/j.sigpro.2014.03.051

CrossRef Full Text | Google Scholar

Hu, J., Wang, Z., Gao, H., and Stergioulas, L. K. (2012). Extended Kalman Filtering with Stochastic Nonlinearities and Multiple Missing Measurements. Automatica 48 (9), 2007–2015. doi:10.1016/j.automatica.2012.03.027

CrossRef Full Text | Google Scholar

Ismail, Z. H., and Jalaludin, N. A. (2016). River Flow and Stage Estimation With Missing Observation Data Using Multi Imputation Particle Filter (MIPF) Method. J. Telecommun. Electron. Comput. Eng. (JTEC) 8 (11), 145–150.

Google Scholar

Jain, S. K., Mani, P., Jain, S. K., Prakash, P., Singh, V. P., Tullos, D., et al. (2018). A Brief Review of Flood Forecasting Techniques and Their Applications. Int. J. River Basin Manage. 16 (3), 329–344. doi:10.1080/15715124.2017.1411920

CrossRef Full Text | Google Scholar

Kang, H. (2013). The Prevention and Handling of the Missing Data. Korean J. Anesthesiol. 64 (5), 402–406. doi:10.4097/kjae.2013.64.5.402

CrossRef Full Text | Google Scholar

Kim, S., Seo, D.-J., Riazi, H., and Shin, C. (2014). Improving Water Quality Forecasting via Data Assimilation – Application of Maximum Likelihood Ensemble Filter to HSPF. J. Hydrol. 519, 2797–2809. doi:10.1016/j.jhydrol.2014.09.051

CrossRef Full Text | Google Scholar

Lee, H.-C., Lin, C.-Y., Lin, C.-H., Hsu, S.-W., and King, C.-T. (2011). A Low-Cost Method for Measuring Surface Currents and Modeling Drifting Objects. IEEE Trans. Instrum. Meas. 60 (3), 980–989. doi:10.1109/tim.2010.2062730

CrossRef Full Text | Google Scholar

Li, D. Z., Wang, W., and Ismail, F. (2014). A Mutated Particle Filter Technique for System State Estimation and Battery Life Prediction. IEEE Trans. Instrum. Meas. 63 (8), 2034–2043. doi:10.1109/tim.2014.2303534

CrossRef Full Text | Google Scholar

Liang, K., Lingfu, K., and Peiliang, W. U. (2012). “Adaptive Gaussian Particle Filter for Nonlinear State Estimation,” in 31st Chinese Control Conference, Hefei, China, July 25-27, 2012, 2146–2150.

Google Scholar

Litrico, X., and Fromion, V. (2009). “Modeling of Open Channel Flow,” in Modeling and Control of Hydrosystems. 1st ed. (London, UK: Springer-Verlag), 17–41.

Google Scholar

Liu, W.-Q., Wang, X.-M., and Deng, Z.-L. (2017). Robust Centralized and Weighted Measurement Fusion Kalman Estimators for Uncertain Multisensor Systems with Linearly Correlated white Noises. Inf. Fusion 35, 11–25. doi:10.1016/j.inffus.2016.08.002

CrossRef Full Text | Google Scholar

Liu, Y., and Gupta, H. V. (2007). Uncertainty in Hydrologic Modeling: Toward an Integrated Data Assimilation Framework. Water Resour. Res. 43, 1–18. doi:10.1029/2006wr005756

PubMed Abstract | CrossRef Full Text | Google Scholar

Mao, J., Ding, D., Song, Y., Liu, Y., and Alsaadi, F. E. (2017). Event-based Recursive Filtering for Time-Delayed Stochastic Nonlinear Systems with Missing Measurements. Signal. Process. 134, 158–165. doi:10.1016/j.sigpro.2016.12.004

CrossRef Full Text | Google Scholar

Mukherjee, A., and Sengupta, A. (2010). Likelihood Function Modeling of Particle Filter in Presence of Non-stationary Non-gaussian Measurement Noise. Signal. Process. 90 (6), 1873–1885. doi:10.1016/j.sigpro.2009.12.005

CrossRef Full Text | Google Scholar

Ogundijo, O. E., Elmas, A., and Wang, X. (2016). Reverse Engineering Gene Regulatory Networks from Measurement with Missing Values. J. Bioinform Sys Biol. 2017 (2017), 2. doi:10.1186/s13637-016-0055-8

CrossRef Full Text | Google Scholar

Pintelon, R., Ugryumova, D., Vandersteen, G., Louarroudi, E., and Lataire, J. (2017). Time-Variant Frequency Response Function Measurement in the Presence of Missing Data. IEEE Trans. Instrum. Meas. 66 (11), 3091–3099. doi:10.1109/tim.2017.2728218

CrossRef Full Text | Google Scholar

Rafieeinasab, A., Seo, D.-J., Lee, H., and Kim, S. (2014). Comparative Evaluation of Maximum Likelihood Ensemble Filter and Ensemble Kalman Filter for Real-Time Assimilation of Streamflow Data into Operational Hydrologic Models. J. Hydrol. 519, 2663–2675. doi:10.1016/j.jhydrol.2014.06.052

CrossRef Full Text | Google Scholar

Rigatos, G. G. (2012). A Derivative-free Kalman Filtering Approach to State Estimation-Based Control of Nonlinear Systems. IEEE Trans. Ind. Electron. 59 (10), 3987–3997. doi:10.1109/tie.2011.2159954

CrossRef Full Text | Google Scholar

Sichangi, A. W., Wang, L., and Hu, Z. (2018). Estimation of River Discharge Solely from Remote-Sensing Derived Data: An Initial Study over the Yangtze River. Remote Sens. 10 (9). doi:10.3390/rs10091385

CrossRef Full Text | Google Scholar

Smith, P. J., Thornhill, G. D., Dance, S. L., Lawless, A. S., Mason, D. C., and Nichols, N. K. (2013). Data Assimilation for State and Parameter Estimation: Application to Morphodynamic Modelling. Q.J.R. Meteorol. Soc. 139 (671), 314–327. doi:10.1002/qj.1944

CrossRef Full Text | Google Scholar

Solonen, A., Hakkarainen, J., Ilin, A., Abbas, M., and Bibov, A. (2014). Estimating Model Error Covariance Matrix Parameters in Extended Kalman Filtering. Nonlin. Process. Geophys. 21 (5), 919–927. doi:10.5194/npg-21-919-2014

CrossRef Full Text | Google Scholar

Tarpanelli, A., Santi, E., Tourian, M. J., Filippucci, P., Amarnath, G., and Brocca, L. (2019). Daily River Discharge Estimates by Merging Satellite Optical Sensors and Radar Altimetry through Artificial Neural Network. IEEE Trans. Geosci. Remote Sensing 57 (1), 329–341. doi:10.1109/tgrs.2018.2854625

CrossRef Full Text | Google Scholar

Tinka, A., Rafiee, M., and Bayen, A. M. (2013). Floating Sensor Networks for River Studies. IEEE Syst. J. 7 (1), 36–49. doi:10.1109/jsyst.2012.2204914

CrossRef Full Text | Google Scholar

Ugryumova, D., Pintelon, R., Vandersteen, G., and Member, S. (2015). Frequency Response Function Estimation in the Presence of Missing Output Data. IEEE Trans. Instrum. Meas. 64 (2), 541–553. doi:10.1109/tim.2014.2342431

CrossRef Full Text | Google Scholar

Wang, S., Fang, H., and Tian, X. (2017). Robust Estimator Design for Networked Uncertain Systems with Imperfect Measurements and Uncertain-Covariance Noises. Neurocomputing 230, 40–47. doi:10.1016/j.neucom.2016.11.035

CrossRef Full Text | Google Scholar

Zhang, X.-P., Khwaja, A. S., Luo, J.-A., Housfater, A. S., and Anpalagan, A. (2015). Multiple Imputations Particle Filters: Convergence and Performance Analyses for Nonlinear State Estimation with Missing Data. IEEE J. Sel. Top. Signal. Process. 9 (8), 1536–1547. doi:10.1109/jstsp.2015.2465360

CrossRef Full Text | Google Scholar

Zhang, X., Khwaja, A. S., Luo, J., Housfater, A. S., and Anpalagan, A. (2014). Convergence Analysis of Multiple Imputations Particle Filters for Dealing with Missing Data in Nonlinear Problems. IEEE J. Sel. Top. Signal. Process. 9 (8), 2567–2570. doi:10.1109/iscas.2014.6865697

CrossRef Full Text | Google Scholar