A High-Precision Method of Stiffness Axes Identification for Axisymmetric Resonator Gyroscopes

Axisymmetric resonators are key elements of Coriolis vibratory gyroscopes (CVGs). The performance of a CVG is closely related to the stiffness and damping symmetry of its resonator. The stiffness symmetry of a resonator can be effectively improved by electrostatic tuning or mechanical trimming, both of which need an accurate knowledge of the azimuth angles of the two stiffness axes of the resonator. Considering that the motion of a non-ideal axisymmetric resonator can be decomposed as two principal oscillations with two different natural frequencies along two orthogonal stiffness axes, this paper introduces a novel high-precision method of stiffness axes identification. The method is based on measurements of the phase difference between the signals detected at two orthogonal sensing electrodes when an axisymmetric resonator is released from all the control forces of the force-to-rebalance mode and from different initial pattern angles. Except for simplicity, our method works with the eight-electrodes configuration, in no need of additional electrodes or detectors. Furthermore, the method is insensitive to the variation of natural frequencies and operates properly in the cases of either large or small frequency splits. The introduced method is tested on a resonator gyroscope, and two stiffness axes azimuth angles are obtained with a resolution better than 0.1°. A comparison of the experimental results and theoretical model simulations confirmed the validity of our method.


Introduction
The family of Coriolis vibratory gyroscopes (CVGs), currently one of the fastest developing types of gyroscope technologies, can be implemented in a wide range of rotation sensing applications, from low-cost consumer markets to high-end military navigations [1][2][3][4][5][6][7]. They can be divided into two broad categories based on the symmetry of its mechanical element: mode-degenerate gyroscopes and mode-nondegenerate gyroscopes [4,6,8]. As a typical kind of mode-degenerate gyroscope, axisymmetric resonator gyroscopes are the most fascinating for their ability to be operated in the whole-angle mode, which provides a fundamentally unlimited dynamical range and a stable scale factor [4,9].
To develop a high-precision CVG, the two most important characteristics of the resonator are a high-quality factor and the high symmetry of the physical parameters. A high-quality factor is a key feature favored by all kinds of resonators, such as optical resonators [10][11][12], nuclear magnetic resonators [13], and mechanical resonators [6]. For all those resonators, a higher-quality factor corresponds to less energy dissipation, and therefore, lower power consumption is needed and less drive force disturbance is included.
An ideal axisymmetric resonator has perfect stiffness symmetry and no damping. As a result of the imperfections of a real axisymmetric resonator, there is both damping and stiffness asymmetry, which are the two main error sources degrading the performance of the gyroscope [8,[14][15][16][17]. Usually, little can be done to reduce the damping asymmetry effectively with repeatability. Stiffness asymmetry leads to a formation of two normal oscillation modes with specific oscillation directions (stiffness axes) and associated natural frequencies. The angle between the two stiffness axes is theoretically determined to be 45 • exactly [1][2][3]. The arbitrary motion of the resonator can be regarded as a superposition of the two principal oscillations [1][2][3]. Stiffness asymmetry can be effectively diminished through both mechanical correction and electrostatic correction [18][19][20][21][22]. To perform either mechanical or electrostatic correction accurately and effectively, it is critical to locate the exact stiffness axes orientations of the resonator.
The conventional approach to obtain the stiffness axes of an axisymmetric resonator is through frequency-sweeping tests or frequency responses analyses [23][24][25][26][27]. Most typically, Ding et al. introduced two direct and simple methods by studying the critical points in the magnitude-frequency responses of the two-dimensional transfer function and by analyzing the peak and the valley values of the beat signals at the readout electrodes, respectively, to obtain the stiffness axes based on the frequency-sweeping and the ringdown approaches [28]. Closkey et al. proposed an experimental system identification method to estimate the orientation of the stiffness axes [29]. They applied modern system identification algorithms to determine the multi-variable input-output model based on the I/O data from the drive electrodes to the sense electrodes. From the model parameter matrices, the stiffness axes are estimated. Schein et al. proposed a parametric modeling technique to estimate the stiffness axes for degenerate modes from the response data [30]. In their approach, they overcame the difficulties confronted by frequency-response-based methods when the frequency split is too small. The modal frequencies, time constants, and stiffness axes can be determined simultaneously in their approach.
This paper provides theoretical analyses and detailed approaches to identify the stiffness axes based upon a measurement of phase differences in signals detected from two orthogonal sensing electrodes for different pattern angles. Firstly, in our approach, the theory is direct and simple; therefore, no time-consuming calculation is needed. Then, our approach shows a weak dependence on the resonant frequency, which means the method works well even when the natural frequency of the resonator has some significant variation due to environmental change. Moreover, our method works well for both large and small frequency splits. Most importantly, our approach focuses mainly on the characterization of the azimuth angle of the stiffness axes. The influences from the modal frequencies, damping, and damping asymmetries are mostly decoupled. Lastly, all external forces are switched off in the measurements; therefore, the errors from inaccurate forces or phases will be minimized.
The theoretical basis of our method is introduced in Section 2, followed by a detailed description of the procedures step by step. Section 3 validates our approach with a resonator gyroscope by comparing the experimental results with the simulations of the theoretical model. Section 4 concludes the paper with a summary of the results.

Materials and Methods
A typical axisymmetric resonator, operating in the n = 2 mode of vibration, can be modeled as a two-dimensional harmonic oscillator. For the ease of illustration, we neglect damping mismatch, assume zero rotation input, and focus on effects of stiffness asymmetry on the dynamics of the system. The equations of motion can be expressed as [8,14] x + 2 where x and y are displacements of the resonator detected by two orthogonal sensing electrodes. As depicted in Figure 1, the angle between two neighboring electrodes in our study is 45 • . The center frequency of the resonator is ω = , with ω L denoting the eigen frequency of the low-frequency axis and ω H the eigen frequency of the high-frequency axis. The frequency split is defined by ∆ω = The angle between the x signal and the high-frequency axis is denoted by θ ω , which is the focus of the paper, as shown in Figure 1. The nominal energy decay time constant is denoted by τ. In this work, we ignore the damping asymmetry of the resonator, because the effect of damping on the dynamics of the resonator is trivial in our method when external control forces are turned off. The equivalent mass m is determined by the structure and material of the resonator.
The external forces f x and f y are forces acting on x and y directions, respectively. The external forces are used for the energy compensation, the quadrature error elimination, and the pattern angle control. To focus on the stiffness asymmetry of the system, we neglect the damping of the system and switch off all external forces temporarily. The equations of motion in this case can be reduced toẍ Furthermore, with the following orthogonal transformation the equations can be easily simplified tö In the new representation, the two modes of vibration are fully decoupled. x and y describe two principal oscillations along two stiffness axes with natural frequencies of ω L and ω H . From Equation (4), the motion of the resonator can be considered as a superposition of the two principal oscillations. The solutions of the two equations can be easily obtained as where the two oscillation amplitudes A L and A H and the two phases φ L and φ H are all completely determined by initial conditions of the resonator. From another perspective, we can set the values of A L , A H , φ L , and φ H by controlling the resonator in different initial states.
Even though the analysis of x and y is simpler and more direct, the signals that can be directly extracted are from the two sensing electrodes x and y. The signals x and y can be calculated by the inverse transformation of (3), from x and y .
Obviously, x and y are completely determined by directions of stiffness axes (θ ω ) and initial conditions (A L , A H , φ L , φ H ) of the resonator. It is possible that by varying initial conditions and looking into the corresponding detected signals, the directions of stiffness axes can be obtained. The core idea in finding the stiffness axes is based on these facts: • If the motion of the resonator has only one principal oscillation mode with nonzero amplitude, which means either A L or A H is zero, then the pattern angle will coincide with one of the two stiffness axes. As a result, the phase variations of the two detected signals will be synchronous, i.e., the two detected signals will have the same phase or the opposite phase, depending on the angles between the pattern angle and the directions of the x and y sensing electrodes. • If both principal vibration modes have nonzero amplitudes, the pattern angle will be between the two stiffness axes. Then, the phase variations of the two detected signals will be asynchronous, i.e., the phase difference between x and y will vary with time.
To facilitate the phase comparison between the two detected signals, we firstly control the resonator in the force-to-rebalance (FTR) mode [31,32]: the two principal vibration modes are forced to oscillate at the same frequency and the phase difference between them is 0, i.e., φ L − φ H = 0. In the FTR mode, the pattern angle can be controlled at an arbitrary direction by setting amplitudes A L and A H to some certain value. Secondly, we switch off all external forces, the detected signals will evolve as expressed in (6). By measuring the phase difference between x and y, we can tell whether the pattern vibration direction coincides with one of the two stiffness axes. The following paragraphs will give a more detailed mathematical description of our approach.
As shown in Figure 2, setting the resonator pattern angle to β, and defining angle α = β − θ ω , the motion of the resonator can be decomposed to two principal vibrations. The amplitudes of the two principal oscillations will be The constant A 0 , a measure of amplitude or energy of the whole pattern vibration, have little influence on the measurement results. However, it is worth noting that A 0 should not be too large so that the system is in the linear region [33]. At time t = 0, we release the FTR control by switching off all external forces; as a result, the two principal oscillations recover their natural frequencies. The two detected signals then will evolve as S x = A 0 cos 2α cos 2θ ω cos ω H t − A 0 sin 2α sin 2θ ω cos ω L t S y = A 0 cos 2α sin 2θ ω cos ω H t + A 0 sin 2α cos 2θ ω cos ω L t (8)

The Amplitudes of Detected Signals versus Time
Firstly, we study the amplitudes of the two detected signals as functions of time. The two detected signals (8) can be easily simplified to with instantaneous amplitudes and phases defined by A x = A 0 2 2 + 2 cos 4α cos 4θ ω − 2 sin 4α sin 4θ ω cos dωt A y = A 0 2 2 − 2 cos 4α cos 4θ ω + 2 sin 4θ ω sin 4α cos dωt ϕ = arctan − sin 2α sin 2θ ω sin dωt cos 2α cos 2θ ω − sin 2α sin 2θ ω cos dωt ϕ = arctan sin 2α cos 2θ ω sin dωt cos 2α sin 2θ ω + sin 2α cos 2θ ω cos dωt Without loss of generality, we set θ ω = 30 • , β = 26 • , f = 7000 Hz, d f = 4 mHz, and give a numerical simulation of S x , A x , S y , and A y in Figure 3. Like most cases, both A L = 0 and A H = 0, the two detected signals S x and S y are additions of two waves with different frequencies, resulting in beats, as shown in Figure 3. From expression (10) and simulation results in Figure 3, it is obvious that the beat frequency is exactly the same to the frequency split d f = dω 2π . The variations of A x and A y versus time are in opposite directions, i.e., when A x reaches the maximum, A y reaches the minimum, and vice versa.
A more detailed simulation of A x and A y with different pattern angles for two cases of θ ω , θ ω = 30 • and θ ω = 70 • , are plotted In Figure 4. In both cases, the variation ranges of A x and A y are closely related to α: the closer the pattern angle is to the stiffness axes, the smaller variation range A x and A y will be; when the pattern angle coincides with either of the stiffness axes, A x and A y will be constant. In theory, we can use this phenomenon to search for the stiffness axes of the resonator. However, the resolution of this approach by observing signal amplitudes variation range is limited in experiments. When the pattern angle is close to the stiffness axes, the amplitude variation range is so small that it is difficult to tell whether the amplitude variation range is increasing or decreasing with small variations of pattern angles.  Experimentally, the amplitudes A x can be obtained easily from detected signal S x by means of conventional phase-locked loop (PLL), as shown in Figure 5.

The Demodulated Signal versus Time
The main information we utilized to identify the azimuth angles of the two stiffness axes is the phase difference between the detected signals from x and y channels, which can be analyzed through demodulation, as described in Figure 5. Using the PLL, we trace the phase and amplitude of the detected signal S x , and then the phase is shifted for π 2 . This shifted signal is then multiplied by the other detected signal S y , and then passed through a low-pass filter. Theoretically, the final demodulated signal, denoted by S yQ , can be calculated as When the resonator is controlled in the FTR mode, the two principal vibration modes are forced to oscillate at the same frequency, i.e., dω = 0, and the demodulated signal S yQ will remain zero. When all the control forces are switched off, the demodulated signal S yQ will evolve as expressed by Equation (11). It is worth noting that Equation (11) shows no direct dependence on θ ω or β but only on their difference α. The slope of S yQ over time is used to tell whether the current β coincides with either stiffness axeṡ The dependence ofṠ yQ on β, in the range of −45 • to 45 • , is plotted in Figure 6 with θ ω = −15 • . At two critical pattern angles when β 0 = θ ω and β 1 = θ ω + 45 • , we havė S yQ = 0. At the high-frequency axis when β = β 0 , the slope ofṠ yQ over β is positive. At the low-frequency axis β = β 1 , however, the slope ofṠ yQ over β is negative. In short, when the pattern angle coincides with either stiffness axes,Ṡ yQ will be zero, and we can distinguish whether the stiffness axis is the low-frequency or the high-frequency axis by observing the slopes ofṠ yQ over β at twoṠ yQ = 0 angles. A detailed procedure for the identification of stiffness axes can be summarized as follows: First, the resonator is controlled in FTR mode with pattern angles β varying from −45 • to 45 • with a step of appropriate angles. At different βs, after the FTR control is stabilized, all control forces are switched off and then the development of the demodulation quantity S yQ is observed for a few seconds; two specific βs, whereṠ yQ = 0, correspond to the two azimuth angles of the two stiffness axes. We can tell either β corresponds to the low-frequency or the high-frequency stiffness axis by measuring slopes ofṠ yQ over β, for which past studies are not discussed in detail.
Even though the damping of the resonator does not have evident influence on the measurement results in this study, it exists in real measurement. For the purpose of comparing with experimental results, an exponential decay of amplitudes is included in the simulations in the following studies.

Results and Discussion
To verify the theoretical analyses and confirm the validity of our method, experiments are conducted for comparison and verification. The experimental setup is displayed in Figure 7a, including a prototype of a hemispherical resonator gyroscope (HRG) and a circuits system. The HRG prototype is composed of a fused silica hemispherical resonator with a high-quality factor and an electrode plate with eight discrete electrodes, as shown in Figure 7b. The electrical connection and control strategy of the circuits system are depicted in Figure 8. The applications of the eight electrodes are clearly shown. The electrodes separated by 180 degrees are connected directly. The x sensing direction is defined as the 0 • direction and the y sensing direction as 45 • . The electrodes with directions orthogonal to the sensing electrodes directions are used for driving. The circuits system consists of a buffer board mounted on the HRG, a mixed-signal board with two analog-to-digital converters (ADC) and two digital-to-analog converters (DAC), and a digital control board, mounted on the top of the mixed-signal board, with both FPGA and ARM functions. The buffer board converts the detected current signals to the appropriate voltage signals, then the mixed-signal board converts the analog voltages to digital signals with two ADCs and generates drive forces with two DACs. The control board implements the forceto-rebalance (FTR) control and the measurement algorithms. The FPGA part performs high-speed demodulation, modulation and DDS processes, and the ARM part performs the other low-speed but more complicated FTR or measurements functions.  Utilizing the proposed procedure, we measured the natural frequencies and the stiffness axes of the resonator, as listed in Table 1. The resolution of the measured stiffness axes angles is better than 0.1 • , which means if we change the β by 0.1 • around each stiffness axis, we can see the opposite sign ofṠ yQ . It should be noticed that the angle between the two stiffness axes obtained is β 1 − β 0 = 43.3 • < 45 • , which is not consistent with the theoretical expectation. This discrepancy could probably be attributed to the fact that the two sensing electrodes are not exactly orthogonal to each other: the gains of the two sensing channels are not the same, or the angle between the two sensing electrodes is not exactly 45 • . This leads to a small deviation in all the measured angles. In this work, we neglect this effect for the moment, and our future work will focus on this problem. Table 1. Experimental results and measured parameters of the HRG; β 0 and β 1 represent directions of high-frequency stiffness axis and low-frequency stiffness axis, respectively. To check the validity of our method, we plotted the amplitudes A x and demodulated signals S yQ as functions of time around the stiffness axes directions based on the experimental measurements and compared the experimental results with the theoretical predictions. For a clear illustration, we use the same parameters as listed in Table 1 for the theoretical simulations, except for a change of β 1 to β 1 = β 0 + 45 • instead of β 1 = β 0 + 43.3 • so as to be consistent with the theory.

Amplitudes versus Time
The amplitude behavior as the function of time of the resonator is presented in Figure 9.

Demodulated Signals versus Time
The demodulated signal S yQ as the function of time, from both the experimental measurements ((b) and (d)) and theoretical predictions ((a) and (c)), for different pattern angles, are shown in Figure 10. The different colors represent different βs. Subfigures (a) and (b) present the results for the pattern angles near the high-frequency stiffness axis, β 0 = −17.8 • . When β < β 0 , we see a negative slope of the demodulated signal S yQ dependence on time in the first few seconds. When β > β 0 , the sign of the slopė S yQ switched from negative to positive. Subfigures (c) and (d) show the behavior of the demodulated signal around the low-frequency stiffness axis. In contrast to the case with the high-frequency stiffness axis, we see a positive slope of S yQ at pattern angles smaller than β 1 and a negative slope at pattern angles larger than β 1 . Again, we see the simulation results fit well with the experiments. These results confirmed the validity of our method from a different perspective.

Conclusions
In conclusion, this paper introduced a novel high-precision identification method of stiffness axes for axisymmetric resonators. Ignoring both the effects of damping asymmetry and the external rotation, the whole motion of the resonator can be decomposed as a superposition of two independent principal oscillations with two different natural frequencies along two stiffness axes. This method is based on comparing the phase differences in the two detected signals. The method was applied on a resonator gyroscope with an electrode plate of eight independent electrodes, and the stiffness axes angles measured show a resolution better than 0.1 • . Finally, a comparison of the detected signal amplitudes A x and demodulated signals S yQ for different pattern angles between the experimental results and theoretical prediction were carried out. A good agreement between them confirmed the validity of the proposed method.
Lastly, even though we take electrostatic sensing and actuation in our experiments on an HRG prototype in this paper, our method is not restricted to these types of actuation or sensing mechanisms and is not limited to HRGs. The core idea behind our method is that the resonator can be excited in any selected pattern angles and can be detected from two non-parallel or non-identical directions. When the resonator is excited in a certain pattern angle, and no clear phase difference deviation is observed in a short time after all the forces are released, the certain angle corresponds to one of the stiffness axes angles.

Data Availability Statement:
The data that support the findings of this study are available from the corresponding author upon reasonable request.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: CVG Coriolis vibratory gyroscopes FTR Force to rebalance PLL Phase-locked loop