Eur. Phys. J. D (2015) 69: 261 DOI: 10.1140/epjd/e2015-60372-4

THE EUROPEAN PHYSICAL JOURNAL D

Regular Article

High precision estimation of inertial rotation via the extended Kalman filter Lijun Liu1,a , Bo Qi2 , Shuming Cheng2,3 , and Zairong Xi2 1 2 3

College of Mathematics & Computer Science, Shanxi Normal University, Linfen 041000, P.R. China Key Laboratory of Systems and Control, Academy of Mathematics and Systems Science, Chinese Academy of Sciences, Beijing 100190, P.R. China University of Chinese Academy of Sciences, Beijing 100049, P.R. China Received 24 June 2015 / Received in final form 9 October 2015 c EDP Sciences, Societ` Published online 24 November 2015 –  a Italiana di Fisica, Springer-Verlag 2015 Abstract. Recent developments in technology have enabled atomic gyroscopes to become the most sensitive inertial sensors. Atomic spin gyroscopes essentially output an estimate of the inertial rotation rate to be measured. In this paper, we present a simple yet efficient estimation method, the extended Kalman filter (EKF), for the atomic spin gyroscope. Numerical results show that the EKF method is much more accurate than the steady-state estimation method, which is used in the most sensitive atomic gyroscopes at present. Specifically, the root-mean-squared errors obtained by the EKF method are at least 103 times smaller than those obtained by the steady-state estimation method under the same response time.

1 Introduction As one of the key sensors of orientation, gyroscopes have been of great importance for practical purposes. They have found wide applications in many areas. Examples include inertial navigation systems [1–3], the stabilization of flying vehicles like radio-controlled helicopters or unmanned aerial vehicles [4,5], and so on. A variety of operating principles have been utilized for sensing changes of orientation, such as mechanical gyroscopes [6,7], fibre optic gyroscopes [8–10], and extremely sensitive quantum gyroscopes [11–14]. For the mechanical and fibre optic gyroscopes, it has been very hard to further improve the precision of the rotation measurements. However, recent technological advancements have enabled quantum gyroscopes to surpass them as the most sensitive rotation sensors. For atomic spin gyroscopes, different spin-ensembles are overlapped in a comagnetometer arrangement to sense the changes of rotation [15]. Compared to the superconducting quantum interface devices [16–18], which are the main competitors, atomic spin gyroscopes have the potential to be significantly cheaper to construct and maintain while also exhibiting better sensitivity. In this paper, we focus on atomic spin gyroscopes and investigate the methods of improving the measurement precision of an unknown inertial rotation rate. The whole measurement process of atomic spin gyroscopes is actually a process of parameter estimation [19–21], and may be divided into three stages. First, the involved spin-ensembles are prepared in an initial state a

e-mail: [email protected]

and evolve under the action of the inertial rotation rate. The second stage concerns the choice of an appropriate measurement applied to extract information about the rotation rate to be estimated. For example, a Faraday modulation technique can be used to measure the optical rotation of the probe beam caused by the alkali-metal atoms. The last step is to associate each experimental result with an estimate of the inertial rotation rate through some rule (e.g. an estimation method). Therefore, the choice of the estimation method in the last step will affect the final actual precision of the ultra sensitive measurement. For atomic spin gyroscopes, the conventional estimation method at present is the steady-state estimation (SSE) method [15,22], which employs the dependency relationship between the steady state signal and the input rotation rate. However, this method cannot work until the spin-ensembles have achieved the steady state. In addition, there exist many approximations in the SSE method which will inevitably affect the final actual precision of the measurement. In this paper, we present an estimation method, the extended Kalman filter (EKF) [23–27], for the third stage of the ultra sensitive measurement. The extended Kalman filter has been used widely in engineering to estimate the state of systems and parameters from noisy measurement signals [28,29]. We introduce this method to atomic spin gyroscopes and demonstrate that by adopting the EKF method as the estimation rule, the practical measurement precision can be greatly improved. Specifically, the rootmean-squared errors obtained by the EKF method can be significantly reduced in comparison with those obtained by the SSE method under the same response time.

Page 2 of 6

The remainder of the paper is organized as follows. In Section 2 we first set up the atomic spin gyroscope model (ASG). We present how to give an estimate of the inertial rotation rate by utilizing the SSE method and the EKF method in Sections 3 and 4, respectively. The estimation errors of the two methods are compared by simulations in Section 5. Section 6 concludes the paper.

2 Atomic spin gyroscopes Atomic spin gyroscopes are based on the comagnetometer arrangement, which in general consists of overlapping spin-ensembles as well as quenching gas (e.g., N2 ) in a spherical glass cell with diameter around 2.5 cm [22]. Here, we describe the K-3 He comagnetometer. K atoms are polarized by a circularly polarized pumping light, while 3 He are polarized through spin-exchange collisions with K atoms [30,31]. The spins of the alkali-metal atoms and the noble gas will evolve under the inertial rotation rate Ω. An off-resonant linearly polarized probe beam which passes through the cell perpendicular to the pump laser can monitor the spin motion of the K atoms. The polarization axis of the probe light rotates by an angle that is proportional to the electron polarization of the alkali-metal atom along the probe direction. This extremely small optical rotation angle can be measured (e.g. by a Faraday modulator), and will yield a measurement of the electron polarization along the probe direction. For the simplicity of analysis, we take the electron polarization of the alkali-metal atom along the probe direction as the measurement signal. In this paper, the z-axis and x-axis are defined as the directions of the pump beam and the probe beam, respectively. The y-axis is perpendicular to the x-z plane. A full characterization of the overlapping atomic spin systems requires density matrix theory. However, for the situation where the spin-exchange rate is much faster than the precession in the magnetic field, the density matrix assumes a spin-temperature distribution, and the coupled dynamics of the alkali-metal atoms and the noble gas can be approximated by the following set of Bloch equations [15] ∂Pe γe = Ω × Pe + (B + λMn Pn + L + be ) × Pe ∂t q 1 en e + (Rpu spu + Rse Pn + Rm spr − Rtot Pe ) , q ∂Pn = Ω × Pn + γn (B + λMe Pe + bn ) × Pn ∂t ne n + Rse (Pe − Pn ) − Rsd Pn . (1) Here, Pe is the electron polarization for the alkali-metal atom, and Pn is the nuclear polarization for the noble gas. Ω is the inertial rotation rate to be measured. γe and γn are electron and 3 He nuclear gyromagnetic ratios, respectively. q is the slowing-down factor depending on the nuclear spin and the electron polarization. B describes the environmental magnetic field. λMe Pe (λMn Pn ) is the total effective magnetic field due to

Eur. Phys. J. D (2015) 69: 261

the presence of the alkali-metal atoms (noble gas) for the 3 He spin ensemble (K atoms). L is the effective magnetic field for K spin caused by the light shift from pump and probe beams. L can be set to zero, because the pump beam is tuned to the center of the optical resonance and the probe beam is linearly polarized. be and bn are anomalous fields coupling to the electron spin and nuclear spin, respectively. In this paper, for simplicity, we assume that be = 0, bn = 0. Rpu and Rm are the pumping rates of the pump and probe lasers, while spu and spr are the mean photon spin of the pump and probe beams. en Rse is the alkali-metal noble-gas spin-exchange rate for ne an alkali-metal atom, and Rse is the spin-exchange rate e e en for a noble gas atom. Rtot = Rm + Rpu + Rsd + Rse is the total spin relaxation rate for an alkali-metal atom, e n where Rsd is the electron spin destruction rate. Rsd is the nuclear spin relaxation rate. In the comagnetometer, the magnetic field along the direction of the pump beam can be decomposed as Bz = Bc + δBz , where |δBz |  Bc , and Bc (≈−λMn Pnz ) is a compensation field which zeros the sensitivity of the measured signal on the magnetic field perpendicular to the pump and probe beams. The measurement signal is Pex . Since there are inevitably noises in measured signals, the gyroscope can only output an approximation rather than the true value of the inertial rotation rate Ω according to the measurement signal Pex and the predetermined estimation method. Different estimation methods may result in different practical precisions of the gyroscope. In the following sections, we employ two different methods, the conventional steady state estimation and the extended Kalman filter, to give the estimate and further compare their effects.

3 Steady state estimation Bloch equations (1) are bilinear, and thus it is difficult to give the analytical relationship between the measured signal Pex and the inertial rotation rate Ω to be detected. For small environmental magnetic fields and rotation rates, the spin-ensembles will achieve equilibrium rapidly. Therefore, the conventional steady state estimation (SSE) method [22] is to establish the dependency relationship between the inertial rotation rate Ω and the steady signal, by setting ∂Pe ∂Pn = = 0. ∂t ∂t For these nonlinear equations, it is still hard to obtain an analytical solution. A series of approximations make the solution tractable, retaining main features of the behavior of the coupled dynamics. Firstly, note that the equilibrium polarization of the spin-ensembles is nearly along the z-axis and thus it is a good approximation to assume that the angles of the polarization vectors Pe and Pn with respect to the z-axis are small enough so that the longitudinal components Pez and Pnz are not affected by the spin-coupling fields

Eur. Phys. J. D (2015) 69: 261

Page 3 of 6

(B, L, Ω, be , bn ). Solving the steady state equations, we have Pez =

n ne Rpu szpu (Rsd + Rse ) n e e ), en ne Rsd (Rpu + Rsd + Rse ) + Rse (Rpu + Rsd

Pnz =

ne Rpu szpu Rse n (R e e . en ne Rsd pu + Rsd + Rse ) + Rse (Rpu + Rsd )

Secondly, Bloch equations (1) can be further simplified when Pez and Pnz are known. Making a series of approximations based on the typical experimental values of the parameters appearing in Bloch equations (see Tab. 3.1 of Ref. [22]), we can find that in the limit of zero Bx , By and δBz , the gyroscope signal is given by:  e Ωy Pez γe Rtot ne − Pex = (1 − Cse ) 2 e 2 2  γ eff Rtot + γe Ωz  γe ne  +Ωx − e (1 + 2Cse ) Ωz γeff Rtot  en  n Rse 1 Rtot , − + γn λMn γe Pez γn Pnz where γeff =

γe γn γe −qγn ,

γ P z Rne

z = Ω

Ωz γeff

,

n Rtot

=

ne Rse

+

n Rsd

From equation (2), it is clearly seen that the measured signal Pex depends on the rotation rate Ωy . However, a lot of approximations are made to obtain equation (2), which will limit the practical precision of the gyroscope. In the next section, we consider another estimation method, the extended Kalman filter, for the atomic spin gyroscope.

4 Extended Kalman filter Instead of the SSE method, we present the EKF method for the atomic spin gyroscope. The EKF method has been widely used in the estimation problems of nonlinear systems. In the following, we first present the EKF method with a general model, and then apply it to estimate the inertial rotation rate in Bloch equations (1). Since measurements are generally made at discrete time instances, we focus on the discrete-time systems. Consider the general nonlinear system model [24]

zk = hk (xk ) + vk .

xk,k = xk,k−1 + Lk (zk − hk (xk,k−1 )) xk+1,k = fk (xk,k )

(4a) (4b)

Lk = Σk,k−1 Hk Γk−1

(4c)

Γk = Hk Σk,k−1 Hk + Rk

(4d)

Σk,k = Σk,k−1 − Lk Hk Σk,k−1

(4e)

Σk+1,k = Fk Σk,k Fk + Gk Qk Gk .

(4f)

and

ne = γne Pez Rese . Especially, if Ωx = 0 and Ωz = 0, the Cse n tot measured signal can be simplified as:   γn γe Ωy Pez x ne 1− (2) Pe = − q − Cse . e γn Rtot γe

xk+1 = fk (xk ) + gk (xk )wk ,

where  denotes the transpose. The initial state x0 is a Gaussian random variable and independent of {wk } and {vk }. The EKF method is used to give an estimate xk,k at time k of the system state xk by using the measurement information up to time k. This estimation procedure is actually a real-time recursive data processing procedure. Denote xk,k and xk,k−1 as estimates of xk at time step k given the measurement set {z0 , z1 , · · · , zk } and the measurement set {z0 , z1 , · · · , zk−1 }, respectively. Given the initial value of the estimate x0,−1 = X0 and the covariance Σ0,−1 = P0 ≥ 0, the extended Kalman filter algorithm can be implemented as [24]

(3a) (3b)

Here, xk is the system state at the time step k, while zk is the measurement output. Nonlinear functions fk (·), gk (·) and hk (·) are sufficiently smooth. {wk } and {vk } are mutually independent, zero mean, white Gaussian processes with covariances E[wk wl ] = Qk δkl and E[vk vl ] = Rk δkl ,

k (x) k (x) Here, Fk = ∂f∂x |x=xk,k , Hk = ∂h∂x |x=xk,k−1 , and Gk = gk (xk,k ). It is worth pointing out that Σk,k can be considered as the approximate conditional covariance of the estimate xk,k [24]. For the atomic spin gyroscope, the rotation rate to be measured is Ωy . To utilize the EKF algorithm (4) to estimate the unknown rotation, we consider the system state as x = (Pex , Pey , Pez , Pnx , Pny , Pnz , Ωy ) . The continuous Bloch equations (1) can be discretized by a small sampling time dt to be the form of equations (3a), and the details of the discretization are given in Appendix A. Here, the system noises w(1) , w(2) and w(3) come from the magnetic fields Bx , By and δBz respectively. These magnetic fields cannot be zeroed accurately, and they are considered as noises with Gaussian distributions whose means are 0. The system output z of equations (3b) corresponds to the measured signal Pex . The performance of the measurement is ultimately limited by quantum fluctuations associated with the alkali-metal atoms and the probe photons, which correspond to the spin projection noise and photon shot noise [32–35], respectively. Both of them are inversely proportional to the number of alkali-metal atoms N being probed. Therefore, we assume that the measurement noise is a Gaussian process, i.e., vk ∼ N (0, σ 2 ) with Rk = σ 2 ∝ 1/N . Through the discretized Bloch equations and the measured signal, we can adopt the EKF algorithm (4) to obtain the estimate xk,k (7) and the approximate conditional covariance Σk,k (7, 7) of the inertial rotation rate Ωy . The initial estimate X0 and the covariance P0 can be determined by the prior information. It is proved in Appendix B that the approximate conditional variance Σk,k (7, 7) of the estimate is monotonically decreasing and will converge with the limit depending on the steady state of the system.

Page 4 of 6

Eur. Phys. J. D (2015) 69: 261 −3

0

10

−7

SSE(t=20 s) EKF(t=0.1 s) EKF(t=1 s) EKF(t=10 s) EKF(t=20 s)

2

−0.0001

0 −2 −4

−0.0002

−2

10

−4

10

10 −3

10

−2

10

−1

0

10

1

10

10

x k, k (7)

0 0.00005559999

−0.2 −0.4

0.00005559980 −4

0

10

10

−2

10

−1

10

5

10 0

10

1

10

−4

10

−5

10

−6

10

−7

10

−8

10

−9

10

−10

10

−11

10

−5

10

−10

10

−4

10

−3

Ω (rad/s)

10

−2

10

y



k,k

(7, 7)

10

−3

root−mean−squared errors (rad/s)

Outputs

x 10

−20

10

−4

10

−3

10

−2

10

−1

10

0

10

1

10

time (s)

Fig. 1. The upper panel is the simulated measured signal Pex with respect to time t. The middle and the lower panels are

y (i.e., xk,k (7)) and its approximate condithe EKF estimate Ω tional variance (Σk,k (7, 7)) with respect to time t, respectively.

5 Numerical simulations The performances of the SSE method and the EKF method for atomic spin gyroscopes are compared as follows. Although the theoretical bound of the parameter estimation error of the gyroscope instrument is ultimately limited by quantum noises, from a practical viewpoint, different estimation strategies will result in different actual estimation errors. A merit quantifier that allows the comparison of different estimation strategies is the rootmean-squared error (RMSE), which can be denoted as:  m 1 

i − Ωy )2 , RMSE = (Ω m i=1 y

yi is the estimate of Ωy based on the ith simuwhere Ω lation experiment, and m is the total number of simulation experiments. This metric characterizes the deviation of the estimated and the real value of the unknown parameter, thus the estimation method will be better when the merit quantifier is smaller. In the following numerical simulations, the values of the parameters appearing in Bloch equations (1) are chosen from typical values in the experiment of [22] (see Tab. 3.1 of [22]). We begin with an ideal case, i.e., the environmental magnetic fields Bx , By and δBz can be zeroed accurately so that Qk = 0 for all k. However, the measurement output noise vk has to be taken into account due to the quantum spin projection noise and the photon shot noise. Recall that the variance of the output noise vk should be inversely proportional to the number of alkali-metal atoms N . Here, we choose Rk = 10−14 (see Tab. 3.1 of [22]). In Figure 1, we take Ωy as 5.56 × 10−5 rad/s. The SSE method cannot work until the Bloch equations (1) achieve the steady state. The corresponding steady-state estimate is 5.833717 × 10−5 at t = 20 s. In contrast, the EKF

Fig. 2. Root-mean-squared errors (m = 200) of the SSE method and the EKF method at different times t (s) for different Ωy (varying from 10−5 to 10−2 (rad/s)). It can be seen that the EKF method is much faster and more accurate than the SSE method.

method can make a real-time estimation of Ωy on the basis of the available measurement information. For example, when t = 0.1 s, the approximate conditional variance

y (xk,k (7)) is less than 10−15 , Σk,k (7, 7) of the estimate Ω and the corresponding estimate (approximate conditional mean) of Ωy is 5.555691 × 10−5 rad/s which almost equals to the true value. When t = 20 s, the estimate of Ωy obtained by the EKF method is 5.559991×10−5 rad/s, which is much more accurate than the SSE method. In order to demonstrate that the EKF method is more accurate than the SSE, we compare the root-meansquared errors of the two methods with different inertial rotation rates, i.e. Ωy , varying from 10−5 to 10−2 (rad/s). The time for Bloch equations (1) achieving the steady state depends on the value of Ωy , which is about 8 s ∼ 10 s for Ωy ∈ (1×10−5, 1×10−2). For convenience, we compare the steady-state performances at time t = 20 s. The simulation results are illustrated in Figure 2. It can be seen that the EKF method is much faster and more accurate than the SSE method. Specifically, the root-mean-squared errors obtained by the EKF method at t = 0.1 s are already smaller than those obtained by the SSE method at t = 20 s. Moreover, if under the same response time (t = 20 s), the root-mean-squared errors obtained by the EKF method are at least 103 times smaller than those obtained by the SSE method. Next, we consider the practical case where Bx , By and δBz are not zeroed accurately. Taking into account the actual experimental conditions, we may assume Bx , By , and δBz are less than 0.05 μG [22]. Here, we treat Bx , By , δBz ∼ N (0, σ 2 ) with σ 2 = 2.89 × 10−16 . We still choose Rk = 10−14 . The simulation results are illustrated in Figure 3. It also shows that the EKF method is much faster and more accurate than the SSE method. Note that the precision of the estimation and the response time are two key performance indices of the gyroscope. From the above analysis, we can conclude that the EKF method is superior to the SSE for the atomic spin gyroscope.

Eur. Phys. J. D (2015) 69: 261

Page 5 of 6

−2

The discretized Bloch equations are

root−mean−squared errors (rad/s)

10

−4

10

SSE(t=20 s) EKF(t=0.1 s) EKF(t=1 s) EKF(t=10 s) EKF(t=20 s)

(1)

zk =

(1) xk

where

−6

10

−8

⎛

(7)

fk (xk ) −10 −4

10

(3)

+ vk

⎞ (1) fk (xk ) ⎜ ⎟ .. fk (xk ) = ⎝ ⎠, .

10

10

(2)

xk+1 = fk (xk ) + gk (xk )(wk , wk , wk )

−3

10

Ω (rad/s)

−2

10

y

Fig. 3. Root-mean-squared errors (m = 200) of the SSE method and the EKF method at different times t (s) for different Ωy (varying from 5 × 10−5 to 10−2 (rad/s)) when Bx , By and δBz are not zeroed accurately. It can be seen that the EKF method is still much faster and more accurate than the SSE method.

6 Conclusions We have presented a simple yet efficient extended Kalman filter method for high precision static and quasi-static inertia estimation. The performances are compared with those obtained by the conventional steady-state estimation method which is being used in the most sensitive atomic gyroscopes at present. Numerical simulations demonstrate that by using our method the root-meansquared error and the response time can be greatly reduced in comparison with those obtained by the SSE method. Future work is to make the EKF method applicable in practical gyroscopes, e.g., we should investigate their responses to time varying rotations. Note that for both the SSE and the EKF methods, the whole process is divided into three stages, and the measurement is performed after the spin-ensembles have evolved for a short time under the action of the inertial rotation input. In addition, the measured signal Pex may have to be collected from averaging measurements on many runs of the experiment under the same initial state. It is worth investigating estimation strategies when other measurement strategies (e.g., continuous measurement without separation into stages) are employed.

Appendix A For the atomic spin gyroscope, the Bloch equations (1) can be discretized by small time step dt and we take dt as 10−5 s in the simulations. Consider Bx , By and δBz as system noises w(1) , w(2) and w(3) , respectively, and vk as the measurement output noise. Let the system state be x = (Pex , Pey , Pez , Pnx , Pny , Pnz , Ωy ) . In order to obtain the discretized Bloch equations similar to the general nonlinear model (3), we rewrite the system state (1) (2) (3) (1) (2) (3) at time k · dt as xk  (xk , xk , xk , yk , yk , yk , Ωky ) .

(A.1) (A.2)

⎛

⎞ (1) gk (xk ) ⎜ ⎟ .. gk (xk ) = ⎝ ⎠. . (7)

gk (xk )

Components of fk (xk ) are   γe (1) (1) (2) (3) y λMn yk + Ωk xk fk (xk ) = xk + q   γe γe (3) (2) Bc + λMn yk + Ωz xk − q q  1 (1) x en (1) e Rpu spu + Rse yk − Rtot xk dt, + q    γe (2) (2) (1) (3) λMn yk + Ωx xk fk (xk ) = xk + − q   γe γe (3) (1) Bc + λMn yk + Ωz xk + q q  1 (2) en (2) e Rpu sypu + Rse dt, + yk − Rtot xk q   γe (3) (3) (1) (2) λMn yk + Ωx xk fk (xk ) = xk + q   γe (2) (1) λMn yk + Ωky xk − q  1 (3) en (3) e Rpu szpu + Rse dt, + yk − Rtot xk q   (4) (1) (2) (3) γn λMe xk + Ωky yk fk (xk ) = yk +   (3) (2) − γn Bc + γn λMe xk + Ωz yk    (1) (1) ne n (1) xk − yk − Rsd + Rse yk dt,    (5) (2) (1) (3) fk (xk ) = yk + − γn λMe xk + Ωx yk   (3) (1) + γn Bc + γn λMe xk + Ωz yk    (2) (2) ne n (2) − Rsd yk dt, + Rse xk − yk (6) fk

(7)

(xk ) =



 (1) (2) γn λMe xk + Ωx yk   (2) (1) − γn λMe xk + Ωky yk    (3) (3) ne n (3) xk − yk − Rsd + Rse yk dt, (3) yk

fk (xk ) = Ωky .

+

Page 6 of 6

Eur. Phys. J. D (2015) 69: 261

Components of g(xk ) are   γe (3) γe (2) (1) , gk (xk ) = 0, xk , − xk q q   γe (3) γe (1) (2) , gk (xk ) = − xk , 0, xk q q   γe (2) γe (1) (3) xk , − xk , 0 , gk (xk ) = q q   (4) (3) (2) gk (xk ) = 0, γn yk , −γn yk ,   (5) (3) (1) gk (xk ) = −γn yk , 0, γn yk ,   (6) (2) (1) gk (xk ) = γn yk , −γn yk , 0 , (7)

gk (xk ) = (0,

0,

0) .

Based on the discretized Bloch equations and the measured signal, the proposed extended Kalman filter can be implemented. The initial estimate X0 and the covariance P0 for the EKF algorithm can be determined by the prior information. In our simulations, we set X0 = [0, 0, 0.5, 0, 0, 0.035, 0] and P0 = [1 × 10−10 , 1 × 10−10 , 1 × 10−14 , 1 × 10−10 , 1 × 10−10 , 1 × 10−14 , 0.1].

Appendix B We retain the notations introduced in Appendix A. By the definitions of Hk , Gk , Fk and equations (A.1)−(A.2), we can obtain the specific expressions of Hk , Gk , and Fk . With the expressions of Fk , Gk and Qk , by equation (4f), we have Σk+1,k (7, 7) = Σk,k (7, 7). Moreover, by equation (4e), we have Σk+1,k (7, 7) = Σk,k−1 (7, 7) −

2 (1, 7) Σk,k−1

Σk,k−1 (1, 1) + Rk

.

Noting that Σk+1,k ≥ 0, it is easy to know that Σk,k (7, 7) is monotonically decreasing and will converge with the limit depending on the steady state of the system. The work was supported by the National Natural Science Foundation of China under Grants Nos. 61134008, 61227902, and 61374092.

References 1. C.W. Tan, S. Park, IEEE Trans. Instrum. Meas. 54, 2520 (2005) 2. J. Bird, D. Arden, IEEE Wirel. Commun. 18, 28 (2011) 3. B. Wang, Z. Deng, C. Liu, Y. Xia, M. Fu, IEEE Trans. Ind. Electron. 61, 2015 (2014)

4. P. Castillo, A. Dzul, R. Lozano, IEEE Trans. Control Syst. Technol. 12, 510 (2004) 5. H.Y. Chao, Y.C. Cao, Y.Q. Chen, Int. J. Control Automat. Syst. 8, 36 (2010) 6. S.E. Alper, T. Akin, J. Microelectromech. Syst. 14, 707 (2005) 7. R. Feng, J. Bahari, J.D. Jones, A.M. Leung, Sens. Act. A 203, 413 (2013) 8. M. Digonnet, S. Blin, H.K. Kim, V. Dangui, G. Kino, Meas. Sci. Technol. 18, 3089 (2007) 9. S.W. Lloyd, V. Dangui, M.J.F. Digonnet, S. Fan, G.S. Kino, Opt. Lett. 35, 121 (2010) 10. C. Clivati, D. Calonico, G.A. Costanzo, A. Mura, M. Pizzocaro, F. Levi, Opt. Lett. 38, 1092 (2013) 11. D.S. Durfee, Y.K. Shaham, M.A. Kasevich, Phys. Rev. Lett. 97, 240801 (2006) 12. J.J. Cooper, D.W. Hallwood, J.A. Dunningham, Phys. Rev. A 81, 043624 (2010) 13. J. Fang, J. Qin, Sensors 12, 6331 (2012) 14. A. Ajoy, P. Cappellaro, Phys. Rev. A 86, 062104 (2012) 15. T.W. Kornack, R.K. Ghosh, M.V. Romalis, Phys. Rev. Lett. 95, 230801 (2005) 16. R.L. Fagaly, Rev. Sci. Instrum. 77, 101101 (2006) 17. B. Abdo, O. Suchoi, E. Segev, O. Shtempluck, M. Blencowe, E. Buks, Europhys. Lett. 85, 68001 (2009) 18. C. Ryu, P.W. Blackburn, A.A. Blinova, M.G. Boshier, Phys. Rev. Lett. 111, 205301 (2013) 19. D. Braun, Eur. Phys. J. D 59, 521 (2010) 20. V. Giovannetti, S. Lloyd, L. Maccone, Nat. Photonics 5, 222 (2011) 21. X. Rong, P. Huang, X. Kong, X. Xu, F. Shi, Y. Wang, J. Du, Europhys. Lett. 95, 60005 (2011) 22. G. Vasilakis, Ph.D. Dissertation, Princeton University (2011) 23. R.E. Kalman, R.S. Bucy, J. Fluids Eng. 83, 95 (1961) 24. B.D.O. Anderson, J.B. Moore, Optimal Filtering (Prentice Hall, New Jersey, 1979) 25. B. Sinopoli, L. Schenato, M. Franceschetti, K. Poolla, M.I. Jordan, S.S. Sastry, IEEE Trans. Automat. Control 49, 1453 (2004) 26. S. Huang, G. Dissanayake, IEEE Trans. Robot. 23, 1036 (2007) 27. C.J. Chiang, J.L. Yang, W.C. Cheng, J. Power Sources 234, 234 (2013) 28. G.L. Plett, J. Power Sources 134, 277 (2004) 29. V.A. Bavdekar, A.P. Deshpande, S.C. Patwardhan, J. Process Contr. 21, 585 (2011) 30. W.C. Chen, T.R. Gentile, T.G. Walker, E. Babcock, Phys. Rev. A 75, 013416 (2007) 31. B. Lancor, T.G. Walker, Phys. Rev. A 83, 065401 (2011) 32. J.L. Sørensen, J. Hald, E.S. Polzik, Phys. Rev. Lett. 80, 3487 (1998) 33. A.A. Clerk, M.H. Devoret, S.M. Girvin, F. Marquardt, R.J. Schoelkopf, Rev. Mod. Phys. 82, 1155 (2010) 34. B. Opanchuk, M. Egorov, S. Hoffmann, A.I. Sidorov, P.D. Drummond, Europhys. Lett. 97, 50003 (2012) 35. M.H. Wimmer, D. Steinmeyer, K. Hammerer, M. Heurs, Phys. Rev. A 89, 053836 (2014)