Estimation and Error Analysis for Optomechanical Inertial Sensors

A sensor model and methodology to estimate the forcing accelerations measured using a novel optomechanical inertial sensor with the inclusion of stochastic bias and measurement noise processes is presented. A Kalman filter for the estimation of instantaneous sensor bias is developed; the outputs from this calibration step are then employed in two different approaches for the estimation of external accelerations applied to the sensor. The performance of the system is demonstrated using simulated measurements and representative values corresponding to a bench-tested 3.76 Hz oscillator. It is shown that the developed methods produce accurate estimates of the bias over a short calibration step. This information enables precise estimates of acceleration over an extended operation period. These results establish the feasibility of reliably precise acceleration estimates using the presented methods in conjunction with state of the art optomechanical sensing technology.


Introduction
Accelerometers are used in countless applications including: satellite, aircraft, robotic, and automotive inertial navigation. Particularly, in the realm of satellite design these devices are used both as navigation equipment [1][2][3][4] and as scientific instrumentation [5][6][7][8][9][10][11]. In both cases the accuracy of the estimated acceleration directly impacts the performance of the system and its ability to meet critical mission requirements. There is presently a vast array of commercially available inertial sensors, developed using a diverse set of physical principles, providing a wide range of precision and sensitivity [12]. Notably, advances in optomechanics have facilitated the development of ultra-sensitive acceleration sensors with a comparatively small form factor [13]. The quality factor, a unitless parameter indicating the sharpness of the resonance peak of the oscillator [14], is of particular interest in quantifying the performance of these sensors. Previous studies indicate that such systems are capable of attaining an exceptionally high quality factor over a range of high and low frequencies, providing a broad dynamic range for acceleration sensing [15][16][17]. The use of highly accurate 3D manufacturing solutions has resulted in sensors with precise operating parameters that are uniquely suited to traditional estimation techniques without reliance on complex error models.
Inertial sensors do not directly observe forcing accelerations or torques, but rather the forces are estimated by assimilating the measurements of the dynamic response of a mechanical system strapped down to the system whose accelerations are to be measured. The fidelity of the sensor model directly impacts the quality of the estimate, and in turn the performance of the navigation system [18][19][20]. Sensor models can never perfectly characterize the behavior of the system and always possess some degree of uncertainty in the established parameters; for this reason these models are intrinsically stochastic. Since the mechanical system's dynamical response to imposed forces is correlated with its own dynamics, which can be non-linear, the observed imposed accelerations are colored by the mechanical system dynamics. Stochastic processes used in the sensor modeling aim at decorrelating the sensor dynamics and associated modeling errors to estimate the imposed acceleration. Characterization of the sensor model's response to these stochastic terms and their impact on performance is an essential component of practical sensor design. This analysis is routinely conducted for angular rate sensors by incorporating external measurements from attitude sensors into a Kalman filter to estimate the gyro biases and their associated steady-state statistics [21][22][23][24]. This analysis has been demonstrated to significantly improve performance in the area of spacecraft attitude determination. Similar procedures are required for estimation of accelerometer biases; namely some source of external information must be provided intermittently to conduct a calibration. This can be done using absolute position estimates from sources, such as GPS [25][26][27], or attitude estimates [28][29][30].
In this paper, a sensor model and estimators for use with state of the art optomechanical inertial sensors [16,17] are presented. Firstly, a method for sensor calibration using a generalized input acceleration is presented using a discrete time Kalman filter to estimate the instantaneous sensor bias. The information from this procedure is then fed forward to the operational phase where estimates of external accelerations are obtained. Two estimation approaches are presented for the task of estimating acceleration, each of which is compatible with the information obtained in the calibration step. These two phases constitute the operational flow of the sensor's estimation program. Results are presented for each of these filters using representative oscillator parameters and simulated measurements in the presence of a sinusoidal external acceleration. It is shown that the calibration filter is able to provide accurate estimates of the instantaneous bias in the presence of an imperfect generalized force over a short window of time. Additionally, the operational filter (for estimates of acceleration) provides precise estimates of external acceleration whose errors follow the bias over extended periods of time. The bias contribution to the acceleration estimate error is well expressed by the associated estimate error variances.

Sensor Model
The sensor is composed of two monolithic fused silica resonators, each of which are made up of a dual-flexure supported proof mass. The dynamics governing the motion of the proof mass undergoing small deflections can be modeled as a perturbed linear harmonic oscillator. The effects of neglected vibrational modes, as well as thermomechanical and optical noise sources are then captured by a Gaussian white noise process and a bias term. The resultant equation of motion is given as: where x is the proof mass's displacement from the equilibrium position, g(t) is the forcing acceleration, ω is the oscillators natural frequency, and ζ is the damping ratio, which is inversely proportional to the quality factor (Q = 1/2ζ). The bias term b(t) is modeled as a Wiener process:ḃ (t) = n u (t) (2) where {n v (t), n u (t)} are uncorrelated, zero-mean Gaussian random processes with autocorrelation functions: where δ(t) is the Dirac delta function. The contributing process noise sources, including: thermomechanical losses and cavity drift have been investigated in references [16,17], respectively. Highly precise manufacturing processes coupled with pressure and temperature control allow for these noise sources to be modeled as Gaussian. The effect of radiation pressure from the optical measurement device is also conservatively introduced to the system using the Wiener process model for the bias. It is worth noting that strictly speaking the contribution of radiation pressure is not Gaussian, however it is commonly linearized for reduced order modeling [31,32]. For this single degree of freedom model the effect of off-axis accelerations are not considered, but could be augmented to a three-axis accelerometer configuration using this model in a straightforward manner. The equations of motion can be expressed in state-space representation as [33]: In this formulation the state is defined as: The continuous time dynamics for the linear time invariant system can be discretized in a straightforward manner [33,34]: where X k = X(t k ), and Φ(t k+1 , t k ) = e A(t k+1 −t k ) is the state transition matrix, which is derived explicitly in Appendix A following the developments of Skelton [35]. It is assumed that g(t) is constant over a sufficiently small interval (t k+1 − t k ). The discrete random process w k encompasses the effects of both n v (t) and n u (t), and is expressed as: The process noise can be sampled using the covariance (Q = E[w k w T k ]) provided in Appendix B. In conjunction with Equation (5), this allows for the state to be determined at discrete time intervals.
Observations of the oscillator are obtained using a highly sensitive laser displacement measuring interferometer which operates at a fixed sampling frequency. These measurements are expressed as a linear function of the state: where H = 1 0 0 is the measurement mapping matrix. The measurement noise (ν k ) made up of shot noise, electric readout noise, and optical frequency noise are taken to be zero-mean Gaussian with variance σ 2 m , and ν k is Gaussian white noise with variance σ 2 m . It is important to note that g(t) is made up of all frequencies of interest, while the observable frequencies of the sensor are limited by the oscillator frequency, the sampling frequency, and the relevant model uncertainties. Therefore, only a component of the forcing acceleration is sensed by the accelerometer.

Calibration Phase: Estimating the Sensor Bias
The accelerometer bias must be monitored to prevent the accumulated build-up of errors in the estimate of acceleration. To do this a calibration step must be conducted on a reasonable basis. One such method is to carefully apply an input to the sensor over a specified time and observe the response of the system. Knowing that the acceleration cannot be perfectly delivered, we define the provided input as: whereḡ(t) is the desired acceleration at time t, and n g (t) is a zero-mean Gaussian process with variance σ 2 g . The discrete sensor dynamics given in Equation (5) can now be augmented with the imperfect input acceleration: where w k is defined in the same way as Equation (7) with the addition of n g (t): Because the input noise n g (t) and the process noise terms {n v (t), n u (t)} are independent, the Q matrix derived in Appendix B would be of the same form, with the sum of σ 2 g and σ 2 v taking the place of σ 2 v . With the dynamics in this form, the discrete-time Kalman filter can be used to sequentially estimate the state using displacement measurements defined in Equation (8) [34,36,37]. The update step, which is conducted when a new measurement is incorporated is given as: where the caret corresponds to the estimate of the indicated quantity and P is the estimate error covariance matrix. After a measurement is incorporated it is propagated over the sampling interval (t k+1 − t k ). The propagation equations are given as: Using this method, an estimate of the sensor bias and the variance of the bias estimate error can be obtained sequentially in the calibration step. Provided a reasonably steady input of an admissible magnitude is applied, the filter will converge to steady-state behavior quite rapidly. With this in mind, the convergent value of P can be found by solving the discrete Riccati equation [38][39][40]: In doing so, a closed-form solution for the expected filter performance can be obtained based on oscillator parameters, sampling frequency, and process noise parameters. This steady-state analysis is important for sensor design and component selection. Accuracies of the displacement-sensing interferometers, process models of the mechanical elements, and the signal processing specifications can be designed using the mathematics and methods of the steady-state analysis.

Operation Phase: Estimating the Forcing Acceleration
The results obtained from the calibration step can now be incorporated to develop procedure for the estimation of the forcing acceleration using a series of displacement measurements defined in Equation (8). Upon completion of the calibration step, it is assumed that the bias estimate (b) is itself unbiased. Under this assertion, the sensor bias is given as: where t 0 is the time of calibration and n b0 is a zero-mean Gaussian random variable with variance σ 2 b0 . Between calibrations the bias estimate is not updated, sob does not change. Because n b0 and n u (t) are both zero mean, the error in the bias estimate between calibrations will remain zero mean; however the variance in the estimate will grow [41]. The statistical quantities governing this growth are given by: This growth in bias uncertainty between calibrations will be shown to degrade the accuracy of the subsequent acceleration estimate as the time since calibration increases. Using this instantaneous definition for the bias term, the discrete dynamics in Equation (5) can be rewritten as: Using this reduced form of the discrete dynamics, two estimators will now be presented.

Least Squares Formulation
The first estimator considered is a moving least squares method using a batch of N + 1 measurements obtained from the laser interferometer, which are defined in Equation (8).
A measurement obtained at n ∈ {0, 1, 2, . . . , N|N ≥ 2} steps ahead of t k is then expressed as the following: where Φ (n) ij and Ψ (n) ij are used as shorthand for Φ ij (t (k+n) , t k ) and Ψ ij (t (k+n) , t k ), respectively. Using this, a set of measurements can then be expressed as a function of the instantaneous oscillator state and forcing acceleration: Now, the minimum variance estimate of x kẋk g k T can be determined from the set of measurements: where the caret corresponds to the estimate of the indicated term and R is the measurement noise covariance. The elements of R are given in indicial notation as: where i ≤ j and δ ij is the Kronecker delta. A closed form solution for the integral in this expression is derived in Appendix B. The measurement noise has an explicit dependence on time, and as a result will increase with the time since last calibration. As a result the certainty in the estimate of acceleration will intuitively degrade during operation, necessitating periodic calibration. Provided the user has some operational requirements on the acceleration estimates, the calibration cycle can be determined explicitly in this manner.

Kalman Filter Formulation
An estimator for the determination of the oscillator state and forcing acceleration can also be developed using a Kalman filter approach. Returning to Equation (17), the dynamics can be expressed as: The new process noise covarianceQ, which has been augmented with the uncertainty in the instantaneous value of the bias, is given by the following: where the integrals contained in the second matrix are derived in Appendix B. The discrete dynamics can then be rearranged to include the instantaneous forcing acceleration as an additional state. In the absence of a model for g(t), the estimate must be driven by the incorporation of new measurements. This can be done by tuning the process noise associated with the channel to an appropriately scaled value (α).
Using these discrete dynamics, the Kalman filter can be run with slight modifications to the update and propagation equations given in Equations (12) and (13), respectively.
Note here that Φ(t k+1 , t k ) provided in Appendix A is identical to matrix in Equation (25), which allows for its use in the propagation equations above. Like the calibration step, the filter is run using displacement measurements of the proof mass given in Equation (8), which each possess measurement noise covariance R = σ 2 m . Figure 1 summarizes the developments for the calibration and operation filter and provides the framework for intended use.

Numerical Simulation Results
The combined calibration and acceleration estimation procedure is implemented and tested using a simulated dataset generated using the sensor model developed in Section 2.1. The oscillator parameters that define the sensor response are provided in Table 1. These parameters correspond to experimental results obtained from benchtesting of a prototype sensor [16,17]. The parameters that define the stochastic processes driving the system {σ v , σ u } are provided in Table 2. Using this information a truth model has been generated, from which the displacement time history has been extracted and corrupted by noise to create a simulated measurement set. These noisy measurements are then provided to the calibration and acceleration estimation programs. Firstly, the calibration step was conducted over a span of 10 min using a constant input acceleration ofḡ = 1 × 10 −6 m/s 2 . The input noise density corresponds to 10% error in delivered acceleration, and is provided in Table 2. The initial estimate of the bias is significantly offset from the true value (1 m/s 2 ). The true and estimated bias during the calibration are plotted with the bias estimate error and corresponding 3σ bounds in Figure 2. The filter is able to converge in the span of a few seconds, producing an unbiased estimate with root mean square error of 1.9816 ng. The converged variance of the bias estimate error (provided in Table 2) is then fed forward to the acceleration estimation procedure in the next step. Table 2. Random process parameters.

Term
Value Units After the calibration phase is implemented, the acceleration estimation step is carried out over an extended operation interval of 8 h. The initial true bias was set to zero, whileb was set as a random drawing from N (0, σ 2 b0 ). The external acceleration was defined as a sinusoid of amplitude 1 × 10 −5 m/s 2 and frequency 1/100 Hz. The results for the error in estimated acceleration and the 3σ bounds are provided in Figure 3 for the moving least squares procedure and in Figure 4 for the Kalman filter. The estimate errors for both filters start as zero mean; however, as time increases the influence of the bias on the estimate error becomes apparent in both cases. This is expected, as the filter is initialized with a constantb for the initial time, and beyond that point the uncertainty in b(t) is reflected by the growth in the covariance. In both cases the estimate error is contained within the plotted 3σ bounds. The acceleration estimates can be deemed reliable until the 3σ bounds become sufficiently large to necessitate calibration, at which point the new bias estimate and σ b0 can be fed forward to continue operation.
To quantify the filter performance in a manner that is independent of the specific realization of the bias random process, the procedure has been conducted for both estimation approaches in a series of 10,000, 1 h trials to generate statistical measures of performance. The results from this Monte Carlo analysis are presented in Table 3.

Discussion
The sensor model developed in Section 2 has been employed to simulate noisy sensor dynamics and measurements. Using these measurements, both the calibration and operation phase estimation procedures have been conducted in series. Figure 2 shows the results obtained in the calibration phase. It is shown that the filter is able to converge to steady-state operation in less than 5 s, in spite of poor a priori knowledge of the bias state. It is shown that the estimate tracks the true value well over the extended 10 min operating window. Figure 3 shows the results from the operation stage using the moving least squares formulation over a period of 8 h without calibration. During this time, the covariance grows significantly as the certainty in the bias estimate provided by the calibration is diminished. The resultant estimate error is shown to be dominated by the effect of the bias, which is plotted in yellow. This same procedure was conducted over a 1 h window for 10,000 Monte Carlo iterations, 100 of which are plotted with the 3σ bounds in Figure 5. The figure shows that the covariance is consistent with the expected results in that the errors are bounded by red curves. An important observation is that the user can determine the required calibration cycle based on required level of accuracy of the acceleration estimates. Figure 4 shows the results from the operation stage using the Kalman filter formulation. The results here mirror those of the least squares formulation. The estimate error is primarily driven by the bias term over the 8 h period, and the estimate error is contained within the 3σ bounds which are plotted in red. The main distinction between the two is the tuning of the parameter α for the process noise covariance. Figure 6 shows the results from 100 Monte Carlo iterations using the Kalman filter. As is shown, the covariance scaling is consistent with the acceleration estimate error results obtained.

Conclusions
Inertial navigation is a crucial component to the future of autonomy across several disciplines. It is required that accurate and precise measurements of acceleration and angular rates be acquired in spite of various sources of error to conduct this procedure effectively. In this paper a stochastic model for novel optomechanical inertial sensors is presented and utilized to derive procedures for the estimation of sensor bias during a calibration step, and forcing accelerations during operation. These procedures are developed to be fully compatible, as outputs from the calibration filter are used directly in the operation phase, thus constituting a full sequence for accurate acceleration measurement. This procedure is demonstrated using representative sensor parameters and simulated process and measurement noise. The simulation results verify that the Kalman filter used for calibration is able to produce a highly accurate estimate of the instantaneous bias in the presence of an imperfect input, and that both operational filters produce precise estimates of acceleration. The acceleration estimation errors induced by the bias as it evolves over the operational window are shown to be fully captured by the covariance; this result is proposed as a quantitative method to determine the required calibration cycle to maintain a desired degree of accuracy. In summary, this work addresses the critical concerns of: how to estimate the forcing acceleration in the presence of a stochastic bias term, how to estimate this bias term during a calibration phase, and how regularly must calibration phases be conducted to retain the desired degree of precision.  Ψ 11 = 1 p 1 ω 2 2p 1 p 2 + e −p 2 dt (p 2 1 − p 2 2 ) sin p 1 (t − t 0 ) − 2p 1 p 2 cos p 1 (t − t 0 ) Ψ 12 = 1 ω 2 1 − 1 p 1 e −p 2 (t−t 0 ) p 1 cos p 1 (t − t 0 ) + p 2 sin p 1 (t − t 0 )

Appendix B. Process Noise Covariance
The process noise vector w k defined in Equation (5) is composed of a Gaussian white noise process, and a random walk. The instantaneous value of this term is given by: For any admissible value of k, w k is a zero-mean Gaussian random vector with covariance Q = E[w k w T k ]. As such, the elements of this covariance matrix are required to define the multivariate distribution from which drawings of w k are taken. By applying the definitions in Equation (3), the elements of Q are given as: