State Estimation Using Dependent Evidence Fusion: Application to Acoustic Resonance-Based Liquid Level Measurement

Estimating the state of a dynamic system via noisy sensor measurement is a common problem in sensor methods and applications. Most state estimation methods assume that measurement noise and state perturbations can be modeled as random variables with known statistical properties. However in some practical applications, engineers can only get the range of noises, instead of the precise statistical distributions. Hence, in the framework of Dempster-Shafer (DS) evidence theory, a novel state estimatation method by fusing dependent evidence generated from state equation, observation equation and the actual observations of the system states considering bounded noises is presented. It can be iteratively implemented to provide state estimation values calculated from fusion results at every time step. Finally, the proposed method is applied to a low-frequency acoustic resonance level gauge to obtain high-accuracy measurement results.


Introduction
Estimating the state of a dynamic system based on noisy sensor measurements is a common problem in sensor methods and applications [1,2]. Mainstream estimation methods all assume that both the system state noise and measurement noise can be modeled as random variables with known statistical properties. The Kalman filter, which supposes both of noises obey Gaussian distributions, is, by far, the most popular method [3]. The basic Kalman filter is only applicable to linear systems. In order to deal with nonlinear cases Bucy and Sunahara proposed the extended Kalman filter (EKF) [4,5]. The EKF uses the first order Taylor expansion technique to linearize state and observation equations, and then obtains state estimations by the Kalman filter. On the other hand, approximation to a state probability distribution of a nonlinear system is, to a great extent, easier and more feasible than a linear approximation to a nonlinear function [6]. Based on this idea, Gordon and Salmond proposed the particle filter (PF) [6]. The performance of the PF is commonly superior to the EKF because it can usually provide more precise information about state posterior probability distribution than does the EKF, especially when it takes a multimodal shape or noise distributions are non-Gaussian [6,7].
The precondition of the above methods is that the noise statistical properties must be known. However, in some practical applications, what engineers can obtain are not precise statistical distributions [8], but ranges of noises. Hence, a group of state estimation methods considering bounded noises, also known as the bounded-error methods, appeared [9][10][11][12]. Assuming that all variables belong to known compact sets, these methods build simple sets, such as ellipsoids or boxes, guaranteed to contain all state vectors consistent with given constraints. For linear systems, some scholars began to 0, A = ∅ (1) Note that the Dempster's combination rule is meaningful only when ∑ B∩C=∅ m 1 (B)m 2 (C) < 1 , i.e., m 1 and m 2 are not totally conflicting. This rule can be used to synthesize uncertain, imprecise or incomplete information coming from different sources.

The Degree of Dependence and the Combination of Dependent Evidence
In DS evidence theory, Dempster's combination rule is the most important tool for computing a new BBA from two BBAs based on two pieces of evidence. This rule requires that the two pieces of evidence must be independent, which is considered to be a very strong constraint and cannot always be met in practice. Wu, Yang and Liu [17] pointed that if there are two pieces of evidence which are partially derived from the same information source, then both of them are mutually dependent. This interpretation concentrates on the connotation of independence conception in evidence combination operation. In this case, Wu, Yang and Liu [17] proposed the energy of evidence concept, and then, deduced the degree of dependence and the dependency coefficient between the two from the energy of the intersection of the two. Based on these notions, the combination of dependent evidence can be realized. Definition 4 (The energy of evidence E). The energy of evidence E, En(E) is defined as: where |A i | is the number of elements in the focal element A i , n(E) is the number of distinct focal elements in E. Obviously, En(E) have some valuable characteristics: (1) if every m(A i ) = 0, namely, m(Θ) = 1, then En(E) = 0 and the evidence E represents no useful information; (2) if |A i | = 1 and m(Θ) = 0, then En(E) = 1 and the E contains the maximum useful information; (3) En(E) ∈ [0, 1]. Suppose that the BBAs of evidence E 1 and E 2 are m 1 and m 2 , respectively, and their focal element sequences are A i and B j . It is possible that some focal elements of E 1 and E 2 are induced by the same information source. In this case, E 1 and E 2 will be dependent, then the energy of the intersection of the two pieces of evidence can be described by: where D ij denotes dependent focal element, |{D ij }| is the number of distinct D ij and the BBA function m is derived from m 1 and m 2 .
The relationship of En(E 1 ), En(E 2 ) and En(E 1 , E 2 ) is illustrated in Figure 1; especially, En(E 1 , E 2 ) = 0 implies the independence between E 1 and E 2 . The value of En(E 1 , E 2 ) measures the dependency of the two pieces of evidence. where Dij denotes dependent focal element, |{Dij}| is the number of distinct Dij and the BBA function m is derived from m1 and m2.
The relationship of En(E1), En(E2) and En(E1, E2) is illustrated in Figure 1; especially, En(E1, E2) = 0 implies the independence between E1 and E2. The value of En(E1, E2) measures the dependency of the two pieces of evidence.  Definition 5 (The degree of dependence between two pieces of evidence). En(E 1 , E 2 ) is defined as the degree of dependence between E 1 and E 2 . Actually, the partial energies En(E 1 ) − En(E 1 , E 2 ) in E 1 and En(E 2 ) − En(E 1 , E 2 ) in E 2 are independent of each other. If energy En(E 1 , E 2 ) is partitioned into two parts, with each part attached to E 1 and E 2 , respectively, as follows: then two corresponding independent pieces of evidence can be generated from E 1 and E 2 . For E 1 , its final independent energy can be calculated as: Similarly: Definition 6 (The dependency coefficient between two pieces of evidence). The dependency coefficient of E 1 to E 2 is defined as: and the dependency coefficient of E 2 to E 1 is defined as: E 1 and E 2 can be modified by R 12 and R 21 , respectively, to obtain their corresponding independent E 1 and E 2 , their BBA functions are given by: Consequently, the requirement of Dempster's rule is met and the combination of E 1 and E 2 can be implemented according to Dempster's rule in (1). Finally, the combination of E 1 and E 2 is indirectly realized by the combination of E 1 and E 2 . Actually, reference [17] gives the decorrelation method to correct E 1 and E 2 by dependency coefficients such that the corrected E 1 and E 2 can be deemed as the independent evidence and combined using Dempster's rule.  [18,19]). A finite support random set on Θ is a pair (F ,m) where F is a finite family of distinct non-empty subsets of Θ and m is a mapping F → [0, 1] and such that ∑ A∈F m(A) = 1.
F is called the support of the random set and m is called a basic belief assignment. Such a random set (F ,m) is equivalent to a mass function in the sense of Shafer.

Extension Principles
The random set (R,ρ) of ζ, which is the image of random relation (F ,m) of ξ through f, is given by extension principles [20][21][22]: where: M is the number of element of F . The summation in Equation (14) accounts for the fact that more than one focal element A i may yield the same image R j through f.
The key of constructing (R,ρ) is to calculate the image of A i through f. If ξ is a continuous variable on Θ, then Θ = R n , F becomes a finite family of distinct non-empty sub-intervals on Θ. In this case, the process of constructing (R,ρ) is given as follows: For each ξ k in ξ, let its marginal random set be (F k ,m k ) and the focal element of (F k ,m k ) be a interval [a k -,a k + ], then the focal element of (F ,m) can be given as: The image of A can be calculated by using the methods of Interval Analysis [19][20][21]; if A is a convex set, then A has 2 n vertices, denoted as v j (j = 1,2, . . . , 2 n ). If function f has certain properties, the Vertex Method can help reduce the calculation time considerably [22]: is continuous in A and also no extreme point exists in this region (including its boundaries), then the value of interval function can be obtained by: Thus, function f has to be evaluated 2 n times for each focal element A. This computational burden can be further reduced if the hypotheses of the following Proposition 2 hold [21].

Proposition 2.
If f is continuous, if its partial derivatives are also continuous and if f is a strictly monotonic function with respect to each ξ k , k = 1, 2, . . . , n, then: There is a case in point.
]. Assume f and its partial derivatives are all continuous. If f is increasing with respect to ξ 1 and ξ 2 , decreasing with respect to ξ 3 respectively, then f has to be calculated only twice for each focal element A, namely, Totally, 2M evaluations of f are needed to obtain complete (R,ρ). Furthermore, the expectation of (R,ρ) is given by [23]: where R j = [r j -, r j + ], j = 1,2,· · · , N, N is the number of focal element R j .

Dynamic System Model under Bounded Noises
The dynamic systems mode constructed by the state and observation equations is as follows: where the relationship between state x k+1 at time k + 1 and state x k at time k is described as function f. The relationship between observation z k+1 at time k + 1 and state x k+1 at time k + 1 is described as function g. v k and w k are bounded additive state noise variable and observation noise variable respectively, which are independent of each other. These two noises can be approximated to triangle possibility distributions [7], noted as π v and π w , respectively, (the noise distributions are identical at each time step), as shown in Figure 2.
where [v a ,v b ] is the support interval of the state noise, v c is the mode of state noise, similarly:   Figure 3 shows the flow of the proposed recursive algorithm. The following steps will be introduced in detail. Step 1: Construct noise evidence to approximate πv and πw Initially, we construct evidence

Recursive Algorithm of State Estimation Based on Extension Principles and Dependent Evidence Fusion
to approximate the possibility distribution πv of state noise variable vk. For any α Î (0, 1], α cut set of πv is [9]: If there exist α0, α1,, αp−1 which satisfy 0 = α0 < α1 << αp−1 < 1, then their corresponding α-cut sets will satisfy , , where p is a positive integer. Take these p α-cut sets as focal elements with nested closed interval forms, then their corresponding BBAs are:  Figure 3 shows the flow of the proposed recursive algorithm. The following steps will be introduced in detail.  Figure 3 shows the flow of the proposed recursive algorithm. The following steps will be introduced in detail. Step 1: Construct noise evidence to approximate πv and πw Initially, we construct evidence

Recursive Algorithm of State Estimation Based on Extension Principles and Dependent Evidence Fusion
to approximate the possibility distribution πv of state noise variable vk. For any α Î (0, 1], α cut set of πv is [9]: If there exist α0, α1,, αp−1 which satisfy 0 = α0 < α1 << αp−1 < 1, then their corresponding α-cut sets will satisfy , , where p is a positive integer. Take these p α-cut sets as focal elements with nested closed interval forms, then their corresponding BBAs are: Step 1: Construct noise evidence to approximate π v and π w . Initially, we construct evidence (F v k , m v k ) to approximate the possibility distribution π v of state noise variable v k . For any α ∈ (0, 1], α cut set of π v is [9]: If there exist α 0 , α 1 ,· · · , α p−1 which satisfy 0 = α 0 < α 1 <· · · < α p−1 < 1, then their corresponding α-cut sets will satisfy where p is a positive integer. Take these p α-cut sets as focal elements with nested closed interval forms, then their corresponding BBAs are: Figure 4 gives an example that when p = 3, and m can be constructed by uniformly cutting α three times. Distinctly, m corresponds to a possibility distribution π that approximates π v . Certainly, a better approximation of the continuous possibility distribution can be obtained by increasing the number p of cut sets, at the expense of higher computational complexity.  Figure 4 gives an example that when p = 3, and m can be constructed by uniformly cutting α three times. Distinctly, m corresponds to a possibility distribution π that approximates πv. Certainly, a better approximation of the continuous possibility distribution can be obtained by increasing the number p of cut sets, at the expense of higher computational complexity. Figure 4. The possibility distribution of state noise and its evidence construction.
It is worth noticing that m is constructed on the condition that all values outside the support interval [va,vb] are completely impossible. However, in practice, the bounds va and vb are commonly given based on available measurement knowledge or real data, so they may not be precise and the values outside [va,vb] may appear. To account for the imprecision of the support interval [7], by discounting m with a small discount rate v, in which, v k m is defined as [9]: In the course of implementing the proposed algorithm, Θ can be replaced by the closed interval [va¢,vb¢], here va¢ >> va and vb¢ >> vb such that the following interval operations can be done easily.
In the same way, we can construct evidence ( , ) using the possibility distribution πw of observation noise variable wk.
Step 2: Obtain state prediction evidence Considering the influence of noise to the state, we construct the state evidence ( , ) x by adding noise to |k k x : as the inputs of state equation by mapping from the inputs to the outputs based on the extension principles in Equations (13) and (14).
Step 3: Obtain observation prediction evidence Taking the state prediction evidence Step 2 as the input of the observation equation based on the extension principles in (13) and (14).
Step 4: Obtain fusion evidence using the possibility distribution πw of wk: After getting observation zk+1 at time k + 1, considering the influence of noise to the observation, we construct the evidence of zk+1 through adding noise to zk+1: It is worth noticing that m is constructed on the condition that all values outside the support interval [v a ,v b ] are completely impossible. However, in practice, the bounds v a and v b are commonly given based on available measurement knowledge or real data, so they may not be precise and the values outside [v a ,v b ] may appear. To account for the imprecision of the support interval [7], constructs (F v k , m v k ) by discounting m with a small discount rate ε v , in which, m v k is defined as [9]: In the course of implementing the proposed algorithm, Θ can be replaced by the closed interval [v a ,v b ], here v a >> v a and v b >> v b such that the following interval operations can be done easily.
In the same way, we can construct evidence (F w k , m w k ) using the possibility distribution π w of observation noise variable w k .
Step 2: Obtain state prediction evidence. E x k+1|k = (R x k+1|k , ρ x k+1|k ) at time k + 1 from state equation. Suppose the estimation result at time k isx k|k . When k = 1,x 1|1 is initialized as real observation z 1 . Considering the influence of noise to the state, we construct the state evidence (F x k , m x k ) ofx k|k by adding noise tox k|k : , ρ x k+1|k ) by mapping from the inputs to the outputs based on the extension principles in Equations (13) and (14).
Taking the state prediction evidence (R x k+1|k , ρ x k+1|k ) in Step 2 as the input of the observation equation g(x k+1 ), we can get E z k+1|k = (F z k+1|k , m z k+1|k ) based on the extension principles in (13) and (14).
Step 4: Obtain fusion evidence.Ê z k+1 = (F z k+1 ,m z k+1|k ) at time k + 1 in observation domain. Firstly, in Step 1, we get evidence (F w k , m w k ) using the possibility distribution π w of w k : After getting observation z k+1 at time k + 1, considering the influence of noise to the observation, we construct the evidence (F z k+1 , m z k+1 ) of z k+1 through adding noise to z k+1 : Secondly, using Dempster s combination rule, we fuse (F z k+1 , m z k+1 ) and (F z k+1|k , m z k+1|k ) to get the fusion evidenceÊ z k+1 = (F z k+1 ,m z k+1 ) in observation domain at time k + 1. As for the relationship between (F z k+1|k , m z k+1|k ) and (F z k+1 , m z k+1 ). The former is obtained by propagatingx k|k from state ); the latter is constructed by adding noise π w (w) to z k+1 . It can be seen that the former completely comes from the state informationx k|k at past time step k which does not use the observation noise w k+1 (π w (w)), but uses the state noise v k (π v (v)). Because w k+1 (π w (w)) and v k (π v (v)) are independent of each other, so it is believed that the former and the latter are also independent of each other. Hence both of them can be directly fused using Dempster s combination rule.
, ρ x k+1|k ) are certainly mutually dependent. Therefore, the combination of dependent evidence must be used for fusing both of them. For the focal elements of (R x k+1 ,ρ x k+1 ) and (R x k+1|k , ρ x k+1|k ) are the closed intervals on real numbers, here we extend the combination of dependent evidence in the discrete frame of discernment introduced in Section 2.2 to that in the continuous frame of discernment (see the corresponding proposition and example in Appendix A). (R x k+1 ,ρ x k+1 ) and (R x k+1|k , ρ x k+1|k ) can be fused using the extended combination of dependent evidence to get state estimation evidence (F x k+1|k+1 ,m x k+1|k+1 ) at time k + 1.
Finally, Pignistic expectation of (F x k+1|k+1 ,m x k+1|k+1 ) is calculated as state estimation valuê x k+1|k+1 by Equation (20). Using state estimation at time k + 1 to do next iteration, we can estimate state at every time step.
In conclusion, as shown in Figure 3, the whole recursive algorithm is actualized under the framework of DS evidence theory. The corresponding evidence in state and observation domains are not only propagated and transformed by the extension principle, but also fused by the Dempster s combination rule and the proposed combination rule for dependent evidence. Especially, fusion procedure can make that the masses focus to those interval focal elements that contain the system state, so as to get the accurate estimation results, which is the main difference from Nassreddine's method under the framework of the interval analysis. In next section, our approach will be applied to liquid level estimation using an industrial level apparatus to show its better performance than possible with Nassreddine's method.

Application to Liquid Level Measurement
Level measurement methods based on sound reflection phenomena have been successfully applied in some areas of process industry (chemical, waste water treatment, petroleum, etc.) because the level is the main monitored process variable used in industrial alarm systems. Ultrasonic measurement methods, with good directivity, convenient operation and so on, have become some of the most commonly used techniques [24]. Their measuring principle is to emit an ultrasound toward a liquid surface and receive the echos, then to calculate the distance from the surface to the acoustic receiver device by multiplying the sound velocity by the round-trip time [25]. However, this method is susceptible to the quality of the instrument itself and environmental noise, which will deteriorate the measurement accuracy. Besides, if the ultrasound encounters foams, residues, deposits, etc., in the measurement process, it is also prone to parasitic reflection, thereby the ultrasound propagation path is changed, which seriously affects the measurement accuracy [26].
On the contrary, low-frequency sound waves have longer wavelength and it is easy to generate the diffraction phenomenon which can effectively overcome the problem of parasitic reflection due to foams, residues, deposits, etc. When a speaker emits sound waves with a uniform change from a frequency f L to a higher frequency f H toward the surface and a microphone receives the corresponding echoes, the generated standing wave signals extracted in the oscilloscope can be used to calculate the height of the liquid level. Kumperščak and Završnik [25,26] used this idea to measure liquid levels. However, they both directly used observations to calculate the liquid level. In practice, if the measurements obtained by using a speaker and a microphone are not precise enough and if the effect of environmental noise is inevitable, then the deviation of the final measurement results will be unacceptable, which is the most common shortcoming in the present level measurement methods.
In our earlier work [27], we have used the Evidential Reasoning(ER) rule to deal with liquid level estimation with bounded noises, but the ER-based method only provides an initial idea for state estimation under the framework of DS evidence theory and only gives precise estimated results when the level length is less than 1.6 m. In order to improve the evidence fusion-based state estimate method, this paper introduces a new information source, Dempster combination rule and evidence dependence conceptions. We construct the state equation and observation equation based on the principle of level measurement using acoustic standing waves, and then use the proposed algorithm to estimate the frequencies of the standing waves, which can be translated into the liquid level height (0 m-10 m). Compared with the direct measurement method and Nassreddine's method, the estimation results verify that our algorithm has obvious advantages and improves the level estimation accuracy.

Acoustic Standing Wave Level Gauge
The structure of an acoustic standing wave level gauge is shown in Figure 5, and mainly consists of a waveguide (a tube), a speaker, a microphone, a thermometer and a controller. When sound waves in the frequency range [f L ,f H ] generated by a signal generator (audio card and speaker) vertically propagate to a surface and echoes appear, superposition of both waves will generate standing waves. Here, y 1 denotes the sound wave generated by speaker and y 2 denotes the echo reflected by the surface: The synthesis wave of y 1 and y 2 can be expressed as: where A is the amplitude of sound wave, P is the frequency of sound wave and L is the distance from the top of the tube to the surface of liquid as shown in Figure 5. From Equation (29), we know that when L and λ have the following relation, the amplitude of synthesis wave reaches the maximum: In this case, this synthesis wave is defined as the standing wave and its wavelength is: where λ k is the wavelength of kth standing wave, f k is the frequency of kth standing wave (kth resonance frequency) in [f L ,f H ]. c is sound velocity, and T is temperature. where λk is the wavelength of kth standing wave, fk is the frequency of kth standing wave (kth resonance frequency) in [fL,fH]. c is sound velocity, and T is temperature. Substituting Equation (30) into Equation (31), we obtain: where, nk is given as [28]: and: Theoretically, in Equation (33), fk+1 − fk = fF, fF is the fundamental resonance frequency and fk = nk fF, nk Î  + (the set of all positive integers) [26,28]. For example, if L = 9.6 m, T = 23.9 °C, and n = 1, then the fundamental resonance frequency can be calculated by (32)

System Model
Firstly, we consider the resonance frequency as the estimated state and construct the corresponding state equation. If we can continuously collect the resonance frequency f k+1 , then we have the following equations: From Equations (32) and (36), obviously, we can get: Consequentially, we can establish the recursive linear state equation and observation equation, respectively: where x k = f k , z k is the observation of f k , w k and v k are independent noise sequences coming from speaker and microphone, respectively, satisfying the conditions: The intervals [v a ,v b ] and [w a ,w b ] denote the boundaries of the state noise and observation noise, respectively. The state noise v k and observation noise w k can be expressed by possibility distributions π v and π w with the support intervals [v a ,v b ] and [w a ,w b ], respectively.
It should be noted that, in theory, n k in (39) should be taken as a positive integer. However, in practice, it can be only calculated by observations z k and z k+1 according to Equation (33). Because of observation imprecision, the calculated n k is commonly not a positive integer, so it should be approximated as: where " • " denotes the operator that "round numbers to the nearest integer".

Liquid Level Estimation Tests
In order to construct the level gauge in Figure 5, we use a low-precision microphone and speaker to emit and receive cosine sound waves, respectively, an electronic thermometer to collect temperature and a PVC tube with a diameter of d = 75 mm to transmit sounds. The estimated level L is the distance from the surface to the speaker platform. The controller transmits sine or cosine waves to drive the speaker to emit the signals vertically to the liquid surface. We use the software AUDIOSCSI (Brothers Studio, Shenzhen, China) which is based on an audio controller (82801HBM-ICH8M with sampling rate 44,100 Hz, Intel Corporation, Santa Clara, CA, USA) to generate sound waves. The frequencies of wave change with uniform speed from f L = 1000 Hz to f H = 2400 Hz in 5 s. Thus, the microphone receives the synthesis waves and sends them to the controller as shown in Figure 6 (L = 4.6 m). It can be seen that there are the frequencies of 39 adjacent standing waves collected by microphone in [1000 Hz, 2500 Hz]. Figure 7 shows the resonance frequency f k (k = 1,2,· · · , 39) extracted from the spectrum of the synthesis waves by the fast detecting algorithm in [28]. In this experiment, set liquid level distance L = 4.6 m, the ambient temperature is 26. For the state noise v k , we use a high-precision oscillograph (TPS2024, Tektronix, Shanghai, China) to receive the cosine sound waves emitted by the audio controller and speaker in range [1000 Hz, 2500 Hz] and calculate errors about 100 frequency points uniformly selected from 1000 Hz to 2500 Hz. As in [7], the bounds of v k , are taken to be plus or minus three times the standard deviation of errors. So the possibility distributions π v of v k can be constructed as in Figure 8. Where the expectation of π v is 0, standard deviation σ v is 0.1, then the support intervals [v a ,v b ] = [−0.3, 0.3], mode v c = 0. Set α 0 = 0, α 1 = 1/3, α 2 = 2/3, we can get three nested closed intervals and their BBAs to approximate π v as in Figure 8.    Table 1.      Table 1.      Table 1.    Table 1.  In the same way, the possibility distributions π w of w k and the corresponding closed intervals and their BBAs can be constructed as in Figure 9, where, σ w = 1.23, w c = −6.9, so [w a ,w b ] = [−10.59, −3.21]. Furthermore, discounting m at rate ε v = 0.05 and approximating as the closed interval [w c − 100σ w ,w c + 100σ w ] = [−129.7, 115.7], we can construct the evidence (F w k , m w k ) of w k shown in Table 2.    , our recursive algorithm presented in Section 3.2 can be used to estimate resonance frequencies at each step k. Figure 10a gives the estimation results of our method and Nassreddine's method, together with the true values and observations (zk). Figure 10b gives the absolute errors between true values and the estimated values of our method, the estimated values of Nassreddine's method, and zk respectively. It can be seen that the estimation accuracy and convergence of our method are better than those of Nassreddine's method because of the focusing function of the proposed fusion procedure for dependent evidence. Finally, we can calculate the estimate level L by Equation (32) according to the estimated resonance frequencies of our method, Nassreddine's method, and zk respectively as shown in Figures  11a,b gives the corresponding absolute values of length estimation errors. Obviously, the more accurate the estimations of resonance frequencies are, the more accurate the estimations of level L are. As our method always provides more accurate estimations of resonance frequencies, it is therefore always superior.    Figure 7, it can be seen that the first observation value of resonance frequency z 1 = 1023.3 Hz. According to Step 2) in Section 3.2, the first estimation resultx 1|1 is initialized as the real observation z 1 . After obtaining (F v 1 , m v 1 ) and (F w 1 , m w 1 ), our recursive algorithm presented in Section 3.2 can be used to estimate resonance frequencies at each step k. Figure 10a gives the estimation results of our method and Nassreddine's method, together with the true values and observations (z k ). Figure 10b gives the absolute errors between true values and the estimated values of our method, the estimated values of Nassreddine's method, and z k respectively. It can be seen that the estimation accuracy and convergence of our method are better than those of Nassreddine's method because of the focusing function of the proposed fusion procedure for dependent evidence.    , our recursive algorithm presented in Section 3.2 can be used to estimate resonance frequencies at each step k. Figure 10a gives the estimation results of our method and Nassreddine's method, together with the true values and observations (zk). Figure 10b gives the absolute errors between true values and the estimated values of our method, the estimated values of Nassreddine's method, and zk respectively. It can be seen that the estimation accuracy and convergence of our method are better than those of Nassreddine's method because of the focusing function of the proposed fusion procedure for dependent evidence. Finally, we can calculate the estimate level L by Equation (32) according to the estimated resonance frequencies of our method, Nassreddine's method, and zk respectively as shown in Figures  11a,b gives the corresponding absolute values of length estimation errors. Obviously, the more accurate the estimations of resonance frequencies are, the more accurate the estimations of level L are. As our method always provides more accurate estimations of resonance frequencies, it is therefore always superior.  Finally, we can calculate the estimate level L by Equation (32) according to the estimated resonance frequencies of our method, Nassreddine's method, and z k respectively as shown in Figure 11a,b gives the corresponding absolute values of length estimation errors. Obviously, the more accurate the estimations of resonance frequencies are, the more accurate the estimations of level L are. As our method always provides more accurate estimations of resonance frequencies, it is therefore always superior.  Table 3. Here, for every values of L, from above to below, Table 3 gives in order the experiment results of our method, Nassreddine's method, and direct measurement method (namely, substituting zk into (32)). It should be noted that the calculation complexity of our algorithm is relatively high, and meanwhile, with the increase of the measured length of level, the corresponding synthetic wave contains more and more resonance frequency points, so the CPU time will increase, so it needs a longer time (the hardware in this test: CPU E8400, CPU Clock Speed 3.00 GHz, RAM 2 GB). But to the situation that liquid level change relatively slowly, our method is applicable. Certainly, the rapid development of data processing capability of computer hardware will make the complexity less of an issue.

Conclusions
The subject of Nassreddine's method is still interval analysis, it introduces evidence, namely belief-function to elaborate bounded noises. In detail, it gives the improved form of noise bounds (the triangular possibility distribution), here the error interval is extended to the evidence construction, namely many interval focal elements with the corresponding belief assignments. Obviously, the latter has more information than the former. Then, it still uses interval arithmetic and constraint-satisfaction techniques to propagate not only interval focal elements, but also its belief    Table 3. Here, for every values of L, from above to below, Table 3 gives in order the experiment results of our method, Nassreddine's method, and direct measurement method (namely, substituting z k into (32)). It should be noted that the calculation complexity of our algorithm is relatively high, and meanwhile, with the increase of the measured length of level, the corresponding synthetic wave contains more and more resonance frequency points, so the CPU time will increase, so it needs a longer time (the hardware in this test: CPU E8400, CPU Clock Speed 3.00 GHz, RAM 2 GB). But to the situation that liquid level change relatively slowly, our method is applicable. Certainly, the rapid development of data processing capability of computer hardware will make the complexity less of an issue.

Conclusions
The subject of Nassreddine's method is still interval analysis, it introduces evidence, namely belief-function to elaborate bounded noises. In detail, it gives the improved form of noise bounds (the triangular possibility distribution), here the error interval is extended to the evidence construction, namely many interval focal elements with the corresponding belief assignments. Obviously, the latter has more information than the former. Then, it still uses interval arithmetic and constraint-satisfaction techniques to propagate not only interval focal elements, but also its belief assignments, hence its performance is slightly better than that of pure interval propagation. However, the subject of our method is DS evidence theory and random set theory and introduces Dempster combination rule and evidence dependence conceptions. Although we still use Nassreddine's evidence construction technique, the random set description of evidence and extension principle of random set are used to obtain state evidence and observation evidence from the defined information sources and to propagate them in the system equations. The main contribution is to realize the fusion of the propagated evidence, in which the degree of dependence and the combination of dependent evidence are further considered.
As a whole, compared with Nassreddine's method, our method increases the state estimation accuracy. The application in liquid level estimation using an industrial level apparatus shows the efficacy of the proposed method. Certainly, it is worth noting that, in a given application, there are some constraint conditions such as the continuous, monotonic and invertible properties of state and observation equations, and state observability. When they cannot be satisfied, the computational burden will inevitably be increased because of the use of additional complex interval operation algorithms or matrix operation algorithms (for multidimensional states) in [21,22]. Hence fast operation algorithms should be studied in the further. On the other hand, although the proposed procedure of evidence fusion can make the masses focus to those interval focal elements that contain the system state, so as to get the accurate estimation results, how to further evaluate the convergence of fusion using available theories is still a problem worthy studying which will promote the usage of evidence theory in state estimation. We can calculate the degree of dependence between the new fusion evidence and the state estimation evidence by Definition 5 and from Equation (A4), we can get: The new fusion evidence is denoted as E 1 , and the state estimation evidence is denoted as E 2 . From Equations (A5) and (A6), we can get the dependency coefficient between E 1 and E 2 :