Research on Magnetohydrodynamic Angular Rate Sensor Denoising for a Space Laser Stabilization Control System

Abstract

The magnetohydrodynamic angular rate sensor (MHD ARS) is a high-bandwidth, high-accuracy sensor that is increasingly used to measure spacecraft harmonic vibration. However, the amplitude of harmonic vibration is usually on the order of microradian to milliradian, and the induced electric potential signal of MHD ARS is only on the order of nanovolt to microvolt, which is easily disturbed by noise. In this paper, an improved method based on autocorrelation with Complete Ensemble Empirical Mode Decomposition with Adaptive Noise (CEEMDAN) and Wavelet Threshold Denoising (WTD) is proposed to denoise the signal of MHD ARS. Firstly, CEEMDAN is used to decompose noisy signals and obtain intrinsic mode functions (IMFs), and autocorrelation is used to determine the relevant modes where the effective signals are located. Then, the improved threshold and thresholding function are used to denoise the relevant modes. Finally, the denoised signal is obtained by combining the denoised relevant modes. In the experiment, noisy MHD ARS signals were recorded in static and dynamic conditions, and the effects of the proposed method and conventional methods were compared. The results of the Allan variance in the static condition and root-mean-square error in the dynamic condition show that the proposed method can effectively overcome the shortcomings of conventional methods and obtain a better denoising effect.

1. Introduction

Free space laser communication and gravitational wave detection are the focus of current research in space exploration. They must achieve a stable laser beam link between the spacecraft [1,2]. Due to the small divergence angle of the laser beam and the narrow receiving field of view [3], the optical system terminal must ensure the precise pointing and tracking of the laser beam. In laser communication, the requirement for laser beam stabilization accuracy is in a few microradians, and in gravitational wave detection, the requirement is in a few $\mathrm{nrad}/\sqrt{\mathrm{Hz}}$ [4]. Spacecraft harmonic vibrations caused by solar array paddle and reaction wheel are the main source of interference in laser beam stability [5]. Measuring the harmonic vibration of a spacecraft can be used for laser beam stabilization control to achieve precise pointing and tracking [6,7]. Researchers have done a lot of work in improving laser beam pointing and tracking accuracy to improve the stability of the laser beam link in optical systems. For example, Scozzafava J. J. et al., proposed a small laser communication structure mounted in a Magnetohydrodynamic Inertial Reference Unit (MIRU). The MIRU provides fine pointing and stabilization [8], where the magnetohydrodynamic angular rate sensor (MHD ARS) provides the MIRU with the precise platform vibration information required for beam stabilization.

MHD ARS can effectively measure the harmonic vibration in the kilohertz bandwidth [9], and its accuracy can reach the nanoradian level. At present, MHD ARS has been applied to some high-precision optical system terminals, such as the Relay Mirror Experiment [10], the Advanced Land Observing Satellite [11], and the Mars Laser Communication Demonstration program [12]. MHD ARS can be combined with gyro and star tracker to form the attitude measurement system. It is the main technical development direction of spacecraft harmonic vibration measurement instruments. However, the amplitude of the harmonic vibration of a spacecraft is typically in the range of a few hundred to a few microradians; then, the induced electromotive force generated by the MHD ARS is only on the order of nanovolt. The weak signal must be amplified to be detected. Noise such as sensor noise, errors caused by analog-to-digital (AD) reference voltage fluctuations, and electromagnetic interference will also be amplified, which seriously affects the detection accuracy. Pinney C. et al. [13] presented a correlation function to measure the noise floor of MHD ARS and proposed a method to predict the drift error using the relationship between the noise floor and drift error. Iwata T. et al. [14] presented an improved Wiener filter to fuse MHD ARS and a low frequency measurement system to reduce the measurement error and noise. Ji Y. et al. [15] presented an adaptive Kalman filter to reduce the low frequency bandwidth error and noise of MHD ARS. These methods require accurate mathematical models and noise characteristics.

In recent years, an adaptive signal processing method, Empirical Mode Decomposition (EMD) [16,17], has been widely studied and has been applied in attitude measurement signal denoising [18]. Complete Ensemble Empirical Mode Decomposition with Adaptive Noise (CEEMDAN) is an improvement on EMD by adding adaptive white noise and calculating a unique residual signal to obtain the intrinsic mode functions (IMFs). CEEMDAN can decompose the signal into IMFs from low to high frequencies in sequence and reduce the reconstruction error [19], but the result contains some noisy IMFs and ineffective IMFs. Therefore, the denoising ability is affected. Researchers proposed some methods for determining the IMFs where the effective signals are located and removing noise that may be contained in some IMFs. Abdelkader R. et al., presented a method of selecting the relevant modes based on the energy content of each mode and the white noise characteristic [20]. However, logarithmic and integral operations are required in this method. Omitaomu O.A. et al., presented a method of determining the cutoff for separating the modes based on mutual information, and denoising the high-frequency group by Wavelet Threshold Denoising (WTD) [21]. Tian Z.W. and Jia Y.F. presented a method for photoplethysmography denoise using wavelet transform to process the noise dominant components obtained from CEEMADN [22]. Fang J.W. et al., presented a denoising method for surface electromyography signals using WTD for the CEEMDAN reconstructed signals [23]. However, these methods do not consider the possible ineffective modes in low-frequency. The shortcomings of conventional WTD also affect the denoising effect.

In this paper, an improved CEEMDAN-WTD method based on autocorrelation is proposed for the harmonic vibration signal measured by MHD ARS. Autocorrelation function (ACF) is used to determine relevant modes, and improve wavelet threshold. First, the modes are obtained by CEEMDAN. The autocorrelation of each mode is used to determine the relevant modes where the effective signals are located. This step can remove some noise and ineffective modes. Then, the relevant modes are denoised using the improved WTD algorithm. The autocorrelation of each wavelet coefficient sequence is used to adjust the wavelet threshold, which can improve the denoising ability of the WTD algorithm and protect the details of the signal. Meanwhile, a new wavelet thresholding function is proposed to overcome the shortcomings of conventional thresholding functions. Finally, the denoised signal is obtained by combining the denoised relevant modes. In static test, the proposed method effectively reduced the quantization noise, angular random walk, and bias instability. In dynamic test, the proposed method improved the signal-to-noise ratio (SNR), and reduced the root-mean-square error (RMSE). Results in both tests show that the proposed method outperforms the conventional methods. Therefore, the proposed method has good prospects for processing spacecraft harmonic vibration signals measured by MHD ARS.

2. Measurement Principle of MHD ARS

Figure 1a shows the structure of MHD ARS, where the bottom ring is a permanent magnet that generates a magnetic field parallel to the sensitive axis. The upper ring is the channel for the conducting fluid. Due to the small viscosity of the conducting fluid, when the vertical plane of the sensitive axis has an angular rate, the conducting fluid will generate a relative angular rate in the conducting fluid channel [11]. At this time, the conducting fluid will cut the magnetic field and generate an induced electric potential. In general, when measuring harmonic vibration on a spacecraft platform, the induced electric potential is only on the order of nanovolt. The output signal is inevitably affected by noise after amplification.

The vibration of a spacecraft platform reaches the optical load through mechanical transmission. According to the working principle of MHD ARS, it should be fixed on the vertical plane of the vibration axis and ensure that the sensitive axis of MHD ARS is parallel to the vibration axis. Figure 1b shows the installation method of MHD ARS on an optical load. The MHD ARS is installed on the back side of the optical load to ensure that the vibration axis is parallel to the sensitive axis of the MHD ARS. The sensitive axis should be as close to the vibration axis as possible to ensure that the MHD ARS vibrates synchronously with the vibration axis.

3. Improved Algorithm

ACF describes the degree of correlation between the values of a random signal taken at any two different moments. The value of ACF of vibration signal decays slowly with increasing time interval, while the value of ACF of noise signal decays rapidly to zero with increasing time interval [24]. The root mean square (RMS) describes the valid value of a physical quantity, which we use to describe the degree of autocorrelation of the signal at different time intervals. A signal with high dispersion is an effective signal. ACF can be used to distinguish between vibration signal and noise.

Figure 2 shows the overall flow of the algorithm. ${\beta }_1$ and ${\beta }_2$ are parameters that can be adjusted. $R_{RMS}^i$ is the RMS of the ACF for each IMF, and $r_{RMS}^j$ is the RMS of the ACF for each wavelet coefficient sequence. $R_{RMS}^i$ and ${\beta }_1$ are used to determine relevant modes. $r_{RMS}^j$ and ${\beta }_2$ are used to improve the threshold. An improved wavelet thresholding function is proposed.

3.1. Determining Relevant Modes Based on Autocorrelation

CEEMDAN obtains each mode by adding white noise and calculating residual at each decomposition stage. Given a discrete signal $x\left(k\right)$, the steps of CEEMDAN can be summarized as follows [25]:

Step1: Add white Gaussian noise to the original signal as follows:

$$x^m\left(k\right)=x\left(k\right)+\epsilon s^m\left(k\right)$$

(1)

where $\epsilon$ is the noise factor and $s^m\left(k\right)$ is the white noise obeying $N\left(0,1\right)$ distribution. $m=1,2,\cdots ,M$, $M$ is the times of adding noise.

Step2: Decompose each $x^m\left(k\right)$ by EMD. In EMD, the local minimum and maximum of each $x^m\left(k\right)$ are calculated by the cubic spline method. The new mode $IM{F}^m$ is the average of the minimum and maximum values [26]. The first mode of CEEMDAN is obtained by averaging the new modes. The first residual is obtained by subtracting this first mode calculated from the previous signal.

$$IMF1\left(k\right)=\frac{1}{M}{\displaystyle \sum _{m=1}^{M}IM{F}^m}$$

(2)

$$c_1\left(k\right)=x\left(k\right)-IMF1\left(k\right)$$

(3)

where IMF1 is the first mode, and $c_1$ is the first residual.

Step3: Repeat all the above steps for the residual and terminate the algorithm when there are no more than two local minimums or maximums in the residual. IMFs and a residual can be obtained as follows:

$$x\left(k\right)={\displaystyle \sum _{i=1}^{K}IMFi}+c$$

(4)

where $i$ is the number of the modes, and $c$ is the residual.

The ACF of a discrete series $y\left(k\right)$ is calculated as follows [24]:

$$R\left(\tau \right)=\frac{1}{n}{\displaystyle \sum _{k=1}^{n}y\left(k\right)\cdot y\left(k+\tau \right)}$$

(5)

where $R\left(\tau \right)$ is the ACF for $y\left(k\right)$ at time interval $\tau$, $n$ is the length of $y\left(k\right)$.

The effective signal is autocorrelated, and its ACF value does not decay rapidly with increasing time intervals. Therefore, its ACF value is highly discrete with respect to the mean value, which means that the RMS value is larger. The RMS of the ACF for each mode $R_{RMS}^i$ is calculated as follows:

$$R_{RMS}^i=\sqrt{\frac{{\displaystyle \sum _{\tau =0}^N{R}^i{\left(\tau \right)}^2}}{N}}$$

(6)

${\beta }_1$ is the set threshold. According to the principle of ACF, when choosing the value of ${\beta }_1$, the amplitude range of the effective signal should be taken into account, as should the possibility of mode mixing [25], a phenomenon in which the effective signal is decomposed into multiple modes, resulting in small values for some of the relevant modes. When $R_{RMS}^i>{\beta }_1$, the IMF is a relevant mode.

3.2. Improved Wavelet Threshold Based on Autocorrelation

Wavelet transform decomposes the original signal into detail coefficients (cD) and approximate coefficients (cA) [27]. The wavelet decomposition process is shown in Figure 3. In the process, the wavelet coefficients of the signal tend to be stable with the increase in the number of decomposition layers, and the wavelet coefficients of the noise gradually decrease [28].

Wavelet transform is a linear transformation [29]. The wavelet coefficients of the signal have obvious autocorrelation, while the wavelet coefficients of the noise have weak or no autocorrelation. The autocorrelation of the wavelet coefficients is used to adjust the threshold.

$R_{RMS}^j$ is the ACF of wavelet coefficients for each frequency band. The RMS of the ACF for wavelet coefficients in different frequency bands $r_{RMS}^j$ is calculated as follows:

$$r_{RMS}^j=\sqrt{\frac{{\displaystyle \sum _{\tau =0}^W{R}^j{\left(\tau \right)}^2}}{W}}$$

(7)

where $W$ is the number of wavelet coefficients, $j$ is the number of the frequency bands.

The Universal Threshold was proposed by Donoho and Johnstone [30]. The Universal Threshold is given as follows:

$$T=\sigma \sqrt{2\ln N}$$

(8)

where $\sigma$ is the noise standard deviation, $N$ is the data length.

The standard deviation of the noise is calculated as follows [31]:

$$\sigma =\frac{median\left(\left|cD\left(k\right)\right|\right)}{0.6745}$$

(9)

where $median\left(\left|cD\left(k\right)\right|\right)$ is the median of all detail coefficients at each scale [32].

We improved the Universal Threshold by exploiting the ACF of the wavelet coefficients. The improved threshold is given as follows:

$$T_j=\left\{\begin{array}{l}{e}^{-a\cdot r_{RMS}^j}\cdot \sigma \sqrt{2\ln N},r_{RMS}^j\ge {\beta }_2\hfill \\ e^{b/r_{RMS}^j}\cdot \sigma \sqrt{2\ln N},r_{RMS}^j<{\beta }_2\hfill \end{array}\right.$$

(10)

where $T_j$ represents the wavelet threshold of different frequency bands. ${\beta }_2$ is the set threshold. According to the principle of ACF, when choosing the value of ${\beta }_2$, these cases should be considered: the amplitude range of the effective signal; mode mixing causes small values of the relevant modes; and the absolute value of the wavelet coefficients of the decomposed effective signal being larger than that of the noise. $a$ and $b$ can adjust the exponent to prevent it from being too large or too small, which would increase the computational effort.

When $r_{RMS}^j>{\beta }_2$, it indicates that there is an effective signal in this frequency band. $e^{-a\cdot r_{RMS}^j}$ and $e^{b/r_{RMS}^j}$ are monotonically decreasing, monotonically increasing exponential functions, respectively. The value of $e^{-a\cdot r_{RMS}^j}$ decreases rapidly as $r_{RMS}^j$ increases, and the value of $e^{b/r_{RMS}^j}$ increases rapidly as $r_{RMS}^j$ decreases. The above-mentioned properties can be used to change the threshold size adaptively in each frequency band. When the autocorrelation of the signal meets the requirement, $e^{-a\cdot r_{RMS}^j}$ is used to adjust the threshold becomes smaller, so that some weak signals can be retained. Otherwise, $e^{b/r_{RMS}^j}$ is used to adjust the threshold to become larger and remove more noise.

3.3. Improved Thresholding Function

Conventional wavelet thresholding functions have obvious shortcomings. The hard thresholding function is discontinuous at the threshold, and the reconstructed signal will generate additional oscillation. The soft thresholding function compresses the wavelet coefficients larger than the threshold, which can easily cause signal distortion. Both of these conventional thresholds lack adaptability [33]. The following hard and soft thresholding function formulas are given as [34].

$${\omega }_{hard}=\left\{\begin{array}{c}\omega ,\hspace{1em}\left|\omega \right|>T\\ 0,\hspace{1em}\left|\omega \right|<T\end{array}\right.$$

(11)

$${\omega }_{soft}=\left\{\begin{array}{cc}\mathrm{sgn}\left(\omega \right)\left(\left|\omega \right|-T\right),& \hspace{1em}\left|\omega \right|\ge T\hfill \\ 0,& \hspace{1em}\left|\omega \right|<T\hfill \end{array}\right.$$

(12)

In order to overcome the shortcomings of conventional wavelet soft and hard thresholding, we improved the thresholding function as follows:

$${\overline{\omega }}_j^k=\left\{\begin{array}{cc}{\lambda }_j^k{\omega }_j^k+\left(1-{\lambda }_j^k\right)\mathrm{sgn}\left({\omega }_j^k\right)\left(\left|{\omega }_j^k\right|-\frac{T_j}{e^{\sqrt{\left|{\omega }_j^k{}^2-T_j{}^2\right|}}}\right),& \hspace{1em}\left|{\omega }_j^k\right|\ge T_j\hfill \\ 0,& \hspace{1em}\left|{\omega }_j^k\right|<T_j\hfill \end{array}\right.$$

(13)

$${\lambda }_j^k=1-e^{-\frac{\sqrt{\left|{\omega }_j^k{}^2-T_j{}^2\right|}}{\alpha }}$$

(14)

where, ${\omega }_j^k$ represents the wavelet coefficient sequences of different frequency bands, ${\overline{\omega }}_j^k$ represents the adjusted wavelet coefficients. ${\lambda }_j^k$ is a scale factor, which can adjust the change trend of the thresholding function between soft and hard thresholds. $\alpha$ is an adjustment parameter, $\alpha >0$.

3.3.1. Continuity Analysis

${\lambda }_j^k$ ensures that the improved thresholding function is continuous at $\pm T_j$. Therefore, the proposed thresholding function can overcome the disadvantages of the hard thresholding function.

$$\underset{{\omega }_j^k\to T_j}{\lim }{\lambda }_j^k=\underset{{\omega }_j^k\to T_j}{\lim }1-e^{-\frac{\sqrt{\left|{\omega }_j^k{}^2-T_j{}^2\right|}}{\alpha }}=0$$

(15)

$$\begin{array}{ll}\underset{{\omega }_j^k\to T_j{}^{+}}{\lim }{\overline{\omega }}_j^k& =\underset{{\omega }_j^k\to T_j{}^{+}}{\lim }{\lambda }_j^k{\omega }_j^k+\left(1-{\lambda }_j^k\right)\mathrm{sgn}\left({\omega }_j^k\right)\left(\left|{\omega }_j^k\right|-\frac{T_j}{e^{\sqrt{\left|{\omega }_j^k{}^2-T_j{}^2\right|}}}\right)\\ & =\underset{{\omega }_j^k\to T_j{}^{+}}{\lim }{\omega }_j^k-T_j=0\end{array}$$

(16)

$$\begin{array}{ll}\underset{{\omega }_j^k\to T_j{}^{-}}{\lim }{\overline{\omega }}_j^k& =\underset{{\omega }_j^k\to T_j{}^{-}}{\lim }{\lambda }_j^k{\omega }_j^k+\left(1-{\lambda }_j^k\right)\mathrm{sgn}\left({\omega }_j^k\right)\left(\left|{\omega }_j^k\right|-\frac{T_j}{e^{\sqrt{\left|{\omega }_j^k{}^2-T_j{}^2\right|}}}\right)\\ & =\underset{{\omega }_j^k\to T_j{}^{-}}{\lim }{\omega }_j^k-T_j=0\end{array}$$

(17)

3.3.2. Deviation Analysis

When $\alpha \to 0$, ${\lambda }_j^k\to 1$, the improved thresholding function is a hard thresholding function. When $\alpha \to \infty$, ${\lambda }_j^k\to 0$, the deviation of the improved thresholding function from the hard thresholding function is shown as follows:

$$\begin{array}{ll}\underset{{\omega }_j^k\to +\infty }{\lim }{\overline{\omega }}_j^k-{\omega }_j^k& =\underset{{\omega }_j^k\to +\infty }{\lim }\mathrm{sgn}\left({\omega }_j^k\right)\left(\left|{\omega }_j^k\right|-\frac{T_j}{e^{\sqrt{\left|{\omega }_j^k{}^2-T_j{}^2\right|}}}\right)-{\omega }_j^k\\ & =\underset{{\omega }_j^k\to +\infty }{\lim }-\left(\frac{T_j}{e^{\sqrt{\left|{\omega }_j^k{}^2-T_j{}^2\right|}}}\right)=0\end{array}$$

(18)

$$\begin{array}{ll}\underset{{\omega }_j^k\to -\infty }{\lim }{\overline{\omega }}_j^k-{\omega }_j^k& =\underset{{\omega }_j^k\to -\infty }{\lim }\mathrm{sgn}\left({\omega }_j^k\right)\left(\left|{\omega }_j^k\right|-\frac{T_j}{e^{\sqrt{\left|{\omega }_j^k{}^2-T_j{}^2\right|}}}\right)-{\omega }_j^k\\ & =\underset{{\omega }_j^k\to -\infty }{\lim }\left(\frac{T_j}{e^{\sqrt{\left|{\omega }_j^k{}^2-T_j{}^2\right|}}}\right)=0\end{array}$$

(19)

$T_j/e^{\sqrt{\left|{\omega }_j^k{}^2-T_j{}^2\right|}}$ ensures that the deviation of the improved thresholding function from the hard thresholding function is zero at infinity.

By choosing the appropriate $\alpha$, the improved thresholding function can be made to converge quickly to the hard thresholding function. This reduces the reconstruction error caused by compressing wavelet coefficients. $\alpha$ makes the algorithm more flexible and can be applied to different working environments.

Figure 4 shows the denoising characteristics of the proposed thresholding function and the conventional thresholding functions. It can be observed that the improved thresholding function approximates the hard thresholding function more and more rapidly as $\alpha$ decreases. The improved thresholding function overcomes the shortcomings of the conventional thresholding functions. Since the improved thresholding function is continuous at $\pm T_j$ and does not produce additional oscillation like the hard thresholding function, it can enhance the smoothness of the denoised signal and effectively protect the details of the signal. Since it does not compress the wavelet coefficient sequence larger than the threshold as the soft thresholding function, it does not cause signal distortion after reconstruction.

4. Experimental Platform and Result Analysis

In this experiment, we acquired the signals of the MHD ARS measurement system in static and dynamic conditions. The experimental platform is shown in Figure 5. The two-axis turntable generates angular vibration on the azimuth axis. The pitch plane is perpendicular to the azimuth axis. Power supply for motors and MAD ARS. AD conversion module and computer for data acquisition and storage. The resolution of the AD module is 0.9 mV. The algorithms are performed using the computer. We fixed the MHD ARS on the pitch plane of the two-axis turntable and placed it as close as possible to the center point of the pitch plane.

4.1. Static Test Analysis

In this test, after energizing the MHD ARS for 30 min, static data were recorded with a sampling frequency of 100 Hz for 1 h. In CEEMDAN, we set the noise factor of the added gaussian white noise to 0.1. To ensure the denoising effect and save on computational costs, we set the times of adding noise to 100 and the maximum number of iterations to 300. Considering the output of the MHD ARS in the experiment is in the millivolt range and the mode mixing, ${\beta }_1$ and ${\beta }_2$ were both chosen as 0.0005. The mother wavelet function was sym6, and the decomposition scale was 5 layers. We set $a$ to 1000 and $b$ to 0.001. To remove as much noise as possible, $\alpha$ was 0.1. The IMFs are shown in Figure 6.

Considering that conventional CEEMDAN-WTD have been applied to signal denoising [20,21,22,23], we used the CEEMDAN, CEEMDAN-hard thresholding, and CEEMDAN-soft thresholding algorithms to compare the effects of the improved method. Figure 7 shows the $R_{RMS}$ of the IMFs. Figure 8 shows the effects of several methods for denoising, the proposed method significantly reduces the base drift.

The Allan variance can be quantified and analyzed as, for example, the quantization noise (Q), angle random walk (N), bias instability (B), rate random walk (K), rate ramp (R) and other noise sources [35]. We performed the Allan variance analysis on the MHD ARS denoised signal under the static condition. As shown in Figure 9, in the curve of the noise signal, the slope of the oblique part is close to −1/2, and the slope of the straight part is close to zero. Is indicates that the main noise of the measurement system is the angle random walk and bias instability. The proposed method has the lowest curve.

The Allan variance results for random noise are listed in

Table 1. Allan variance analysis results of MHD ARS signal in static condition.

Methods	Q (°)	$\mathbf{N}(°/\sqrt{\mathbf{h}})$	B (°/h)	K (°/h3/2)	R (°/h2)
Noisy MHD ARS signal	7.171 × 10−4	2.088 × 10−3	0.205	3.151	18.056
CEEMDAN	2.225 × 10−3	6.289 × 10−4	0.337	6.539	36.058
CEEMDAN-Soft Thresholding	1.185 × 10−3	4.941 × 10−4	0.242	4.687	25.829
CEEMDAN-Hard Thresholding	2.230 × 10−3	8.403 × 10−4	0.386	7.474	41.206
Proposed Method	5.400 × 10−5	6.466 × 10−5	0.117	2.249	12.370

. As can be seen, the quantization noise is reduced by 10 times as compared to the noisy signal, the angle random walk is reduced by 30 times, and the bias instability is reduced by 1.75 times. The method proposed in this paper has advantages over conventional methods in the static test.

4.2. Dynamic Test Analysis

The working bandwidth of the MHD ARS is 2–1000 Hz. We used the azimuth motor of the turntable to generate angular vibrations of different frequencies to simulate the harmonic vibrations of a spacecraft. In order to verify the effectiveness of the proposed method under the dynamic condition, several single-frequency signals were tested. We acquired the angular vibration signals with different frequencies of 5 Hz, 10 Hz, and 15 Hz with a sampling frequency of 2 kHz and filtered the noisy MHD ARS signals.

The angular vibration can be expressed as follows:

$$\theta \left(t\right)=A\sin \left(2\pi ft\right)$$

(20)

where $A$ is the amplitude of the angular vibration and $f$ is the frequency. The angular rate of the angular vibration can be expressed as follows:

$$w\left(t\right)=2\pi f\cdot A\cos \left(2\pi ft\right)$$

(21)

In this test, the amplitude of the angular vibration is at the microradian level, and the output of the MHD ARS is at the millivolt level. In CEEMDAN, mode mixing results in relatively small magnitudes for the relevant modes. For the above reasons, the parameters were the same as the settings of the static experiment. We used the angular rate measured by an optical encoder system mounted on the azimuth axis as the true angular rate. We acquired single-frequency angular vibration signals at 5 Hz, 10 Hz, and 15 Hz with different angular rates.

As shown in Figure 10, the vibration signals at each frequency contain three amplitudes. At 5 Hz, the angular rates of the vibrations are 0.300°/s, 0.200°/s, and 0.100°/s. At 10 Hz, the angular rates of the vibrations are 0.150°/s, 0.100°/s, and 0.050°/s. At 15 Hz, the angular rates of the vibrations are 0.100°/s, 0.050°/s, and 0.025°/s. Figure 11 shows the IMFs and frequency spectrum of the 0.100°/s signal at 15 Hz.

It is observed that the denoised signal obtained by the proposed method is much smoother. Figure 11 shows the IMFs of the 0.100°/s signal at 15 Hz. Figure 12 shows the $R_{RMS}$ of the relevant modes. As shown in Figure 11, mode mixing appears in the relevant modes, which results in a larger magnitude of IM7 than the other relevant modes. It is due to the fact that noise affects the selection of the local minimum and maximum values, the average values are obtained from the combination of the noise and the real signal [24]. Noise in the relevant modes IMF5 and IMF6 can be observed.

It can be seen from the frequency spectrum in Figure 11 that IMF10 and IMF11 are ineffective modes in low-frequency without the 15 Hz signal, and IMF1-IMF4 are noise modes in high-frequency.

RMSE and SNR are commonly used indicators to evaluate the effect of signal denoising [36].

$$RMSE=\sqrt{\frac{1}{N}{\displaystyle \sum _{k=1}^N{\left(x\left(k\right)-\overline{x}\left(k\right)\right)}^2}}$$

(22)

$$SNR=10\lg \left(\frac{P_S}{P_N}\right)$$

(23)

where $x(k)$ is the true vibration signal of the turntable, $\overline{x}(k)$ is the denoised signal. $P_S$ is the signal power, and $P_N$ is the noise power.

In

Table 2. RMSE and SNR results of MHD ARS signals at 5 Hz.

Methods	Angular Rate (°/s)	RMSE	SNR (dB)	SNR Gain (%)
Noisy MHD ARS signal	0.300	2.564 × 10−3	31.3682	-
	0.200	2.424 × 10−3	28.3353	-
	0.100	1.452 × 10−3	26.8497	-
CEEMDAN	0.300	2.454 × 10−3	31.7504	1.22
	0.200	1.613 × 10−3	31.8698	12.47
	0.100	1.217 × 10−3	28.3814	5.70
CEEMDAN-Soft Thresholding	0.300	2.453 × 10−3	31.7515	1.22
	0.200	1.613 × 10−3	31.8705	12.48
	0.100	1.217 × 10−3	28.3872	5.72
CEEMDAN-Hard Thresholding	0.300	2.382 × 10−3	32.0061	2.03
	0.200	1.613 × 10−3	31.8703	12.48
	0.100	1.120 × 10−3	28.5281	6.25
Proposed Method	0.300	2.005 × 10−3	33.5042	6.81
	0.200	1.333 × 10−3	33.5232	18.31
	0.100	7.462 × 10−4	32.6329	21.54

,

Table 3. RMSE and SNR results of MHD ARS signals at 10 Hz.

Methods	Angular Rate (°/s)	RMSE	SNR (dB)	SNR Gain (%)
Noisy MHD ARS signal	0.150	1.638 × 10−3	29.2376	-
	0.100	1.306 × 10−3	27.6868	-
	0.050	1.309 × 10−3	21.6461	-
CEEMDAN	0.150	1.165 × 10−3	32.2010	10.14
	0.100	9.543 × 10−4	30.2935	9.41
	0.050	9.012 × 10−4	24.8867	14.97
CEEMDAN-Soft Thresholding	0.150	8.007 × 10−4	35.4563	21.27
	0.100	7.434 × 10−4	30.3009	9.44
	0.050	7.697 × 10−4	26.2572	21.30
CEEMDAN-Hard Thresholding	0.150	1.114 × 10−3	32.5879	11.46
	0.100	9.264 × 10−4	30.6676	10.77
	0.050	8.911 × 10−4	24.9846	13.55
Proposed Method	0.150	5.513 × 10−4	38.6975	23.36
	0.100	5.820 × 10−4	34.7052	25.34
	0.050	5.992 × 10−4	28.4324	31.35

and

Table 4. RMSE and SNR results of MHD ARS signals at 15 Hz.

Methods	Angular Rate (°/s)	RMSE	SNR (dB)	SNR Gain (%)
Noisy MHD ARS signal	0.100	1.615 × 10−3	25.8416	-
	0.050	1.717 × 10−3	19.2888	-
	0.025	1.165 × 10−3	16.6340	-
CEEMDAN	0.100	1.419 × 10−3	26.9644	4.34
	0.050	1.140 × 10−3	22.8417	18.42
	0.025	8.549 × 10−4	19.3249	15.64
CEEMDAN-Soft Thresholding	0.100	1.294 × 10−3	27.7629	7.43
	0.050	1.140 × 10−3	22.8421	18.42
	0.025	8.548 × 10−4	19.3252	16.18
CEEMDAN-Hard Thresholding	0.100	1.405 × 10−3	27.0507	4.68
	0.050	1.131 × 10−3	22.9154	18.80
	0.025	8.507 × 10−4	19.3676	16.43
Proposed Method	0.100	1.198 × 10−3	28.4383	10.05
	0.050	8.862 × 10−4	25.0325	29.78
	0.025	6.224 × 10−4	22.0809	32.75

, the RMSE and SNR of vibrations of different angular rates at three frequencies are shown. The RMSE is reduced by 46.58%, and the SNR is improved by 32.75% for the 0.025°/s signal at 15 Hz. The proposed method has the lowest RMSE and the highest percentage of SNR gain for different rates of angular vibration signals at each frequency. It is observed that as the rate decreases, the percentage of SNR gain increases in the improved algorithm. This is because as the rate decreases, the effect of noise on the relevant modes increases. The improved threshold and thresholding function can be used to better advantage.

The harmonic vibration of a spacecraft has high-frequency with low amplitude and low-frequency with high amplitude characteristics [37,38]. To further verify the effectiveness of the proposed method in denoising the MHD ARS dynamic signal, we acquired mixed-frequency angular vibration signals with different amplitudes at 5 Hz, 10 Hz, and 15 Hz.

As shown in Figure 13, the mixed-frequency signal with high amplitude consists of the high amplitude of each frequency. At 5 Hz, the amplitudes of the angular vibrations are 100 $\mathsf{\mu }\mathrm{rad}$, 75 $\mathsf{\mu }\mathrm{rad}$, and 50 $\mathsf{\mu }\mathrm{rad}$, the corresponding angular rates are 0.180°/s, 0.135°/s, and 0.090°/s. At 10 Hz, the amplitudes of the angular vibrations are 25 $\mathsf{\mu }\mathrm{rad}$, 20 $\mathsf{\mu }\mathrm{rad}$, and 15 $\mathsf{\mu }\mathrm{rad}$, the corresponding angular rates are 0.090°/s, 0.072°/s, and 0.054°/s. At 15 Hz, the amplitudes of the angular vibrations are 10 $\mathsf{\mu }\mathrm{rad}$, 8.5 $\mathsf{\mu }\mathrm{rad}$, and 5 $\mathsf{\mu }\mathrm{rad}$, the corresponding angular rates are 0.054°/s, 0.046°/s, and 0.027°/s.

The proposed method has a smoother denoising effect on the mixed-frequency signal with a high amplitude, as shown in Figure 13. Figure 14 and Figure 15 show the IMFs of the mixed-frequency signal with a high amplitude and the $R_{RMS}$ of the relevant modes. The signal is decomposed into multiple modes. Mode mixing occurs in the relevant modes, as shown in Figure 14, which results in a larger magnitude of IM8 than the other relevant modes. Noise in the relevant mode IMF6 can be observed.

It can be seen from the frequency spectrum in Figure 13 that IMF11 is an ineffective mode in low-frequency without the mixed-frequency signal, and IMF1-IMF5 are noise modes in high-frequency.

In

Table 5. RMSE and SNR results of MHD ARS mixed-frequency signals.

Methods	Angular Rate (°/s)	RMSE	SNR (dB)	SNR Gain (%)
Noisy MHD ARS signal	0.180, 0.090, 0.054	1.912 × 10−3	30.7586	-
	0.135, 0.072, 0.046	1.718 × 10−3	29.3816	-
	0.090, 0.054, 0.027	1.306 × 10−3	28.3960	-
CEEMDAN	0.180, 0.090, 0.054	1.740 × 10−3	31.5782	2.66
	0.135, 0.072, 0.046	1.457 × 10−3	30.8117	4.87
	0.090, 0.054, 0.027	1.009 × 10−3	30.6355	7.89
CEEMDAN-Soft Thresholding	0.180, 0.090, 0.054	1.435 × 10−3	33.2485	8.09
	0.135, 0.072, 0.046	1.211 × 10−3	32.4143	10.32
	0.090, 0.054, 0.027	8.126 × 10−4	32.5149	14.51
CEEMDAN-Hard Thresholding	0.180, 0.090, 0.054	1.683 × 10−3	31.8686	3.61
	0.135, 0.072, 0.046	1.418 × 10−3	31.0493	5.68
	0.090, 0.054, 0.027	9.795 × 10−4	30.8916	8.79
Proposed Method	0.180, 0.090, 0.054	1.329 × 10−3	33.9170	10.27
	0.135, 0.072, 0.046	1.091 × 10−3	33.3205	13.41
	0.090, 0.054, 0.027	7.618 × 10−4	33.0754	16.48

, the RMSE and SNR of the MHD ARS mixed-frequency signals with different amplitudes are shown. The RMSE is reduced by 41.65% and the SNR is improved by 16.48% for mixed-frequency signals with a high amplitude. It is observed that the proposed method has the lowest RMSE and the highest percentage of SNR gain for different amplitudes of mixed-frequency angular vibration signals. As with single-frequency signals, the percentage of SNR gain increases as the amplitude decreases in the improved algorithm.

From the dynamic experimental results, it can be seen that the proposed method can effectively improve the SNR and the results are smoother for both single-frequency and mixed-frequency signals.

5. Conclusions

In this paper, an autocorrelation-based CEEMDAN-WTD algorithm is proposed to denoise the static and dynamic signals of the MHD ARS. Firstly, CEEMDAN was used to decompose noisy signal and obtain IMFs, and the autocorrelation of each mode was used to determine the relevant modes. Compared with conventional CEEMDAN, this step could remove ineffective modes in low-frequency and noise IMFs in high-frequency. Then, the autocorrelation of each wavelet coefficient sequence was used to adjust the wavelet threshold. Compared with conventional WTD, this step could remove more noise and retain weak signals. Meanwhile an improved thresholding function was proposed to overcome additional oscillation and reconstruction error in the conventional WTD. Finally, the denoised signal was obtained by combining the denoised relevant modes. Compared with the conventional method, the proposed method could achieve a better denoising effect.

In the static test, the proposed method outperformed the conventional methods. The quantization noise, angle random walk, and bias instability were effectively reduced. In dynamic test, the proposed method improved the SNR and reduced the RMSE more than conventional methods. The denoised signal produced by the proposed method was smoother. As the rate decreased, the effect of noise on the relevant modes increased, the improved threshold and thresholding function achieved a better denoising effect. The tests verified the effectiveness of the proposed method for denoising the angular vibration signal measured by MHD ARS. The proposed method has good prospects in laser communication, gravitational wave detection, and others. Future work is needed to verify the denoising of the angle vibration signals in the three axes and implement it in hardware.