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.

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.

The Bayesian theorem used by the DA method is represented as follows (Crisan and Doucet, 2002). p ( X t | W 1 : t ) = p ( W t | X t ) p ( X t | W 1 : t − 1 ) p ( W t | W 1 : t − 1 ) (14) where X t is the system state at time t, W t is the observation at time t, p ( X t | W 1 : t ) is the posterior probability of state, X at time t given observation, W from time one to time t, p ( W t | X t ) is the likelihood function of state X at time t given observation W at time t, p ( X t | W 1 : t − 1 ) is the prior probability of state X at time t given observation W from time one to time t-1, p ( W t | W 1 : t − 1 ) 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 ( W t | W 1 : t − 1 ) = ∫ p ( W t | X t ) p ( X t | W 1 : t − 1 ) 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 ( W t | X t ) . This function is represented as follow (Mukherjee and Sengupta, 2010) p ( W t | X t ) = p ( W t − g ( X t ) ) (16) where W t and X t are defined in Eq. 14, and g ( X t ) 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 { 0 :     Observation   is   missing   from   sensor   k   at   time   t 1 :     Observation   is   available   from   sensor   k   at   time   t 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 Ξ t , while the available information set Ψ 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) Ξ t j ∼ φ ( Ξ t | Ψ 0 : t ) = ∑ i = 1 N ω ˜ t i p ( Ξ t | X ˜ t i )     for   i = 1 , … , N   and   j = 1 , … , M (17) where Ξ t represent all missing observations at time t, Ψ 0 : t represent all available observation from time 0 to time t, { ω ˜ t i , X ˜ t i } 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 = { Ξ t j , Ψ t } (18) where Ξ t j represent all missing observation during jth imputation and time t, and Ψ t represents all available observation at time t. The imputed sets are used in determining the posterior probability density that represent as p ( X t | Ψ 0 : t ) = ∫ p ( X t | W 0 : t − 1 , Ψ t ) p ( Ξ t | Ψ 0 : t ) d Ξ t (19) where X t is the system state at time t, W 0 : t − 1 is the complete observation, Ξ t and Ψ 0 : t are defined in Eq. 17, and Ψ t is defined in Eq. 18. Considering the Monte Carlo approximation, the probability density can be written as p ( X t | Ψ 0 : t ) ≈ 1 M ∑ j = 1 M p ( X t | W 0 : t − 1 , U t j ) (20) where M is defined in Eq. 17, U t j is defined in Eq. 18, and X t , Ψ 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 ( X t | W 0 : t − 1 , U t j ) ≈ ∑ i = 1 N ω t j , i δ ( X t − X t j , i ) (21) where X t j , i is the system state at ith particle and jth imputation at time instance t, ω 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 ( X t | Ψ 0 : t ) ≈ 1 M ∑ j = 1 M ∑ i = 1 N ω t j , i δ ( X t − X t j , i ) (22) where Ψ 0 : t , M, N are defined in Eq. 17, X t is defined in Eq. 19, and X t j , i , ω 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).

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:

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.