The optical Frequency Scanning Interferometry (FSI)-based absolute ranging technology has broad application prospects in many fields, such as multi-satellite formation, large-scale equipment assembly and manufacturing, and automatic driving^[1-9]. The typical light source of the FSI ranging system is an External Cavity Diode Laser (ECDL), which continuously tunes the optical frequency by changing the length of the external cavity through piezoelectric ceramics or motors. However, due to the nonlinear response of the displacement tuning element, the optical frequency change of the ECDL output is nonlinear with time under the linear drive signal input^[10-12]. The nonlinearity in optical frequency scanning of ECDL leads to the frequency time-varying characteristics of the FSI interference signal, which in turn affects the phase extracting accuracy of the interference signal, and finally leads to the decrease of FSI ranging accuracy^{[2, 11, 13-14]}. Therefore, reducing the influence of optical frequency scanning nonlinearity on interference signals is essential to improve the ranging accuracy of FSI.

The current research methods for the nonlinearity in optical frequency scanning of ECDL can be summarized into the following two categories: one type of method is to realize the early suppression of frequency scanning nonlinearity by monitoring the optical frequency change information and using the feedback loop to control the optical frequency scanning rate online. For example, in 2009, Roos P A etal..^[15] used the double phase-locked loop method to control the optical frequency scanning rate of ECDL online, reducing the influence of scanning nonlinearity on the frequency of interference signals, and achieving a ranging accuracy of 86 nm in the range of 1.5 m. This method suppresses the scanning nonlinearity before the optical path forms interference by introducing a feedback control loop. However, the introduction of feedback control loop increases the system complexity, and the loop bandwidth also limits the increase in optical frequency scanning rate; another type is to eliminate the influence of frequency scanning nonlinearity on the phase extracting accuracy of interference signal by sampling and post-processing the frequency time-varying interference signal. For example, in 2014, SHI G etal..^[2] constructed an all-fiber Mach-Zehnder reference interferometer with a 30 m optical path difference, which uses the generated reference interference signal to resample the measurement interference signal at equal optical frequency intervals, thus greatly reducing the influence of scanning nonlinearity on the interference signal measurement and improving the spatial resolution of the ranging system by about 30 times. However, in such resampling methods, the reference interference signal used for resampling needs to satisfy the Nyquist sampling theorem, which leads to the limitation of the system measurement range by the optical path difference of the reference interferometer. In addition, the introduction of the reference interferometer also increases the complexity of the ranging system.

Both of the above two types of methods increase the system complexity in different extent, and limit the measurement bandwidth and measurement range of the ranging system. To this end, in 2018, DENG W etal..^[16-18] proposed that Empirical Mode Decomposition (EMD) combined with Hilber Transform (HT) to perform phase extracting of nonstationary interference signals, which effectively reduced the influence of the optical frequency scanning nonlinearity on the phase extracting accuracy of interference signals. However, traditional EMDs have component mode aliasing when decomposing and reconstructing nonstationary signals, which limits the further improvement of phase extracting accuracy. Therefore, in this paper, the non-stationary signal mode decomposition reconstruction method based on noise-assisted data analysis method-Complementary Ensemble Empirical Mode Decomposition (CEEMD) is used to decompose and reconstruct the measurement interference signal, improve the modal aliasing caused by the traditional EMD method, and further reduce the influence of the optical frequency scanning nonlinearity.

As shown in Fig. 1, the FSI ranging system outputs a frequency scanning laser from ECDL, which is divided into two beams at the beam splitter (BS1) after being collimated: one beam enters the Fabry-Perot (F-P) Etalons, and the transmitted F-P signal is generated by photodetector (PD1); the other enters the Michelson interference light path and is divided into reference light and measurement light by the beam splitter (BS2). After being reflected by the reflecting prism RR1 and RR2, the interference signal is formed on the photodetector (PD2).

In FSI ranging system, the output optical frequency of ECDL $\upsilon (t)$ can be expressed as:

$ \upsilon (t) = {\upsilon _0} + \beta \cdot t \quad,$ (1)

where ${\upsilon _0}$ is the initial optical frequency output by the laser, and $\beta $ is the optical frequency scanning rate. According to the principle of light wave interference, the interference signal can be expressed as:

$ I(t) = {I_{\rm{r}}} + {I_{\rm{m}}} + 2\sqrt {{I_{\rm{r}}} \cdot {I_{\rm{m}}}} \cos \left[ {2{\text{π}} {\upsilon _0}\tau + 2{\text{π}} \beta \cdot \tau \cdot t} \right]\quad, $ (2)

where $ {I_{\rm{r}}} $ is the light intensity of the reference light, ${I_{\rm{m}}}$ is the light intensity of the measurement light, and $\tau $ is the delay of the measurement light relative to the reference light in the time domain. According to Eq. (2), the instantaneous phase of the interference signal $ \varphi $ is expressed as:

$ \varphi = 2{\text{π}} {\upsilon _0}\tau + 2{\text{π}} \beta \cdot \tau \cdot t \quad.$ (3)

If the optical frequency variation of ECDL is $\Delta \upsilon $ within the period $\Delta t$=t₂-t₁, the corresponding phase difference $\Delta \varphi $ can be expressed as follows:

$ \Delta \varphi = \varphi ({t_2}) - \varphi ({t_1}) = 2{\text{π}} \cdot \tau \cdot \beta \cdot \Delta t \quad.$ (4)

In combination with equation (1), the phase difference $\Delta \varphi $ can be further rewritten as follows:

$ \Delta \varphi = 2{\text{π}} \cdot \tau \cdot \left[ {\upsilon ({t_2}) - \upsilon ({t_1})} \right] = 2{\text{π}} \cdot \tau \cdot \Delta \upsilon\quad. $ (5)

When the distance to be measured is $L$, the delay $\tau $ can be expressed as:

$ \tau = \frac{{2n \cdot L}}{c} \quad,$ (6)

where $n$ is the refractive index of air and $c$ is the velocity of light. The expression of distance L to be measured can be obtained by connect equations (5) with (6):

$ L = \frac{{c \cdot \tau }}{{2n}} = \frac{{c \cdot \Delta \varphi }}{{4{\text{π}} \cdot n \cdot \Delta \upsilon }}\quad. $ (7)

It can be seen from Equation (7) that the ranging accuracy of the system depends on the extraction of the phase difference $\Delta \varphi $ of the interference signal and the calculation of the optical frequency variation range $\Delta \upsilon $. The optical frequency variation range can be obtained by counting the number of F-P signal peaks in Fig. 1, and after selecting the starting and ending F-P peaks, the relationship between the optical frequency variation range $\Delta \upsilon $ and the F-P peak number $r$ is as follows:

$ \Delta \upsilon = (r - 1) \cdot {FSR}\quad, $ (8)

where FSR is the Free Spectrum Range (FSR) of the Fabry Perot etalon, that is, the corresponding optical frequency difference between adjacent peaks. Therefore, solving the interference signal phase difference $\Delta \varphi $ corresponding to the scanning range $\Delta \upsilon $ is the key to determine the ranging accuracy of FSI.

In fact, the interference signal intercepted by the F-P signal includes an integer period signal segment and a fractional period signal segment, and the corresponding phase difference between two adjacent peaks in the integer period signal segment is $2{\text{π}} $, so the integer periodic phase difference can be directly solved by the number of peaks of the interference signal. As shown in Figure 2, in order to improve the solving efficiency of the interference signal phase difference, the interference signal phase difference is divided into three parts, that is, integer phase difference $\Delta {\varphi _{\rm{i}}}$, starting decimal phase difference $\Delta {\varphi _{\rm{s}}}$ and cutoff decimal phase difference $\Delta {\varphi _{\rm{e}}}$. When the phase segment is automatically intercepted, from the moment ${t_{\rm{s}}}$ of the starting F-P peak, the segment is intercepted to the right to the moment ${t_{{\rm{p}}1}}$ of the first wave peak , which corresponds to the starting decimal phase difference interval; similarly, from the moment ${t_{\rm{e}}}$ at which the F-P peak is terminated, the segment is intercepted to the left to the moment ${t_{{\rm{p}}2}}$ at which the last wave peak is, which is the interval of the cut-off decimal phase. The interval from ${t_{{\rm{p}}1}}$ to ${t_{{\rm{p}}2}}$ is the integer phase interval. The integer phase difference $\Delta {\varphi _{\rm{i}}}$ can be expressed as:

$ \Delta {\varphi _{\rm{i}}} = (\delta - 1) \cdot 2{\text{π}}\quad. $ (9)

The integer phase difference is determined by the peak number $\delta $ of the interference signal, the decimal phase difference is $\Delta {\varphi _{\rm{s}}}$, and $\Delta {\varphi _{\rm{e}}}$ is solved by HT:

$ \left\{ \begin{gathered} \Delta {\varphi _{\rm{s}}} = {\varphi _{\rm{u}}}({t_{{\rm{p}}1}} - {t_{\rm{s}}})/T \\ \Delta {\varphi _{\rm{e}}} = {\varphi _{\rm{u}}}({t_{\rm{e}}} - {t_{{\rm{p}}2}})/T \\ \end{gathered} \right. \quad,$ (10)

where ${\varphi _{\rm{u}}}$ is the phase unwrapping function after HT, and $ T $ is the period of the interference signal. Thus, the phase of the interference signal can be expressed as:

$ \Delta \varphi = \Delta {\varphi _{\rm{i}}} + \Delta {\varphi _{\rm{s}}} + \Delta {\varphi _{\rm{e}}} \quad.$ (11)

At this time, the scanning frequency range $\Delta \upsilon $ and the phase difference $\Delta \varphi $ of its corresponding interference signal can be substituted into Eq. (7) to get the distance $L$.

However, the displacement tuning element PZT of ECDL has hysteresis and creep characteristics, resulting in a nonlinear relationship between the optical frequency change output by the ECDL and the time under the linear driving signal input. The optical frequency scanning rate $\;\beta $ is no longer a fixed value, and the interference signal period T becomes a time variable, and an error will occur when calculating the decimal phase using HT. Therefore, this paper proposes to use the CEEMD-HT algorithm to calculate the decimal phase value point-by-point in the intercepted interference signal segment, and then calculate the phase change of the intercepted interference signal, so as to avoid the impact of scanning nonlinearity on its phase calculation.

Firstly, the phase solving method of HT in FSI interference signal is derived, and for any signal $x(t)$, its HT can be expressed as:

$ H[x(t)] = {\mathcal{F}^ - }[\mathcal{F}(x(t))\delta (\omega )]\quad, $ (12)

where $ \mathcal{F}[ \cdot ] $ and ${\mathcal{F}^ - }[ \cdot ]$ are the Fourier transform and the inverse Fourier transform, and $\delta [\omega ]$ is the window function of HT, which can be expressed as:

$ \delta (\omega ) = - j{{\rm{sgn}}} (\omega ) = \left\{ \begin{gathered} - j{\text{ }}(\omega > 0) \\ 0{\text{ (}}\omega = 0{\text{)}} \\ j{\text{ (}}\omega < 0{\text{)}} \\ \end{gathered} \right.\quad. $ (13)

According to Eq. (2), the $ I(t) $ expression of FSI interference signal can be simplified as follows:

$ I(t) = a(t) + b(t)\cos \varphi (t)\quad, $ (14)

where $ a(t) $ is the DC component of the interference signal, $ b(t) $ is the amplitude modulation of the interference signal, and the phase information of the interference signal is hidden in $ b(t)\cos \varphi (t) $. Use a suitable high pass filter to filter out the low-frequency DC term $ a(t) $, and get the term $ b(t)\cos \varphi (t) $ containing the phase information of the interference signal. The interference signal can be expressed as follows:

$ {I_p}(t) = {\mathcal{F}^ - }[\mathcal{F}[I(t)] \cdot W(\omega )] = b(t)\cos \varphi (t)\quad, $ (15)

where $ W(\omega ) $ is the window function of the high pass filter. In addition, BEDROSIAN E^[17]etal.. have made a detailed discussion on the HT in the form of product, and obtained an important conclusion: the HT of the product of low-pass signal and high-pass signal is equivalent to the product of the HT of low-pass signal and high-pass signal. Based on this conclusion, the HT of the interference signal after high-pass filtering can be expressed as:

$ H[{I_p}(t)] = b(t)H[\cos \varphi (t)] = b(t)\sin \varphi (t) \quad.$ (16)

Based on the above analysis, the HT is used to construct the analytic function $ \hat I(t) $ of the interference signal:

$ \hat I(t) = {I_p}(t) + iH[{I_p}(t)]\quad. $ (17)

Combining equations (15) and (16), we can get:

$ \hat I(t) = b(t)\cos \varphi (t) + ib(t)\sin \varphi (t)\quad. $ (18)

The phase function $\varphi (t)$ of the interference signal can be expressed as:

$ \varphi (t) = \arctan \left[ {\frac{{{{\rm{Im}}} (\hat I(t))}}{{{{\rm{Re}}} (\hat I(t))}}} \right] = \arctan \left[ {\frac{{b(t)\sin \varphi (t)}}{{b(t)\cos \varphi (t)}}} \right]\quad. $ (19)

According to Eq. (19), the phase change of interference signal $I(t)$ from time ${t_1}$ to ${t_2}$ can be expressed as:

$ \Delta \varphi = {\varphi _{\rm{u}}}({t_2}) - {\varphi _{\rm{u}}}({t_1})\quad. $ (20)

It should be noted that the Signal-to-Noise Ratio (SNR) of the interference signal is a key factor affecting the accuracy of HT phase solution^[16]. The interference signal is affected by the measurement environment perturbation, photodetector noise and other factors, especially under long-distance measurement conditions, and its SNR is often low, so here the EMD is used to perform adaptive decomposition of the interference signal and reconstruct a high-quality interference signal to improve the interference signal phase extracting accuracy. In the signal decomposition, EMD does not depend on the basis function, and the adaptive decomposition is performed completely according to the extreme point distribution of the signal itself, which decomposes the complex signal containing multiple frequency components into the sum of multiple single frequency components step by step. For the analysis and processing of complex signals with non-stable characteristics, the specific process of EMD decomposition is described as follows:

(1) Firstly extract the upper and lower extreme points of the original signal $s(t)$, and fit the upper and lower extreme points with triple spline fitting to obtain the upper envelope signal ${e_{\rm{u}}}(t)$ and the lower envelope signal ${e_{\rm{d}}}(t)$.

(2) Find the mean envelope signal $m(t)$ of the upper and lower envelopes, remove $m(t)$ from the initial signal $s(t)$, we obtain the signal $c(t)$.

(3) Treat $c(t)$ as a new $s(t)$ signal and repeat steps (1)-(2) until the latest decomposed $c(t)$ satisfies the decomposition termination condition^[17], at which point $c(t)$ is considered a first-order component $IM{F_1}$.

(4) Remove $IM{F_1}$ from $s(t)$, record the remaining signal as $r(t)$, treat $r(t)$ as a new $s(t)$ signal, repeat steps (1)-(3), and obtain the components $IM{F_i}$ of each order in turn.

The decomposition of nonstationary signal $s(t)$ can be expressed as:

$ s(t) = \sum\limits_{i = 1}^n {IM{F_i}} + {r_n}(t)\quad, $ (21)

where $ IM{F_1} $, $IM{F_2},\cdots, $$IM{F_n}$ are the components of each order obtained by EMD, also known as the Intrinsic Modal Function (IMF), and Fig. 3 shows a schematic diagram of the decomposition reconstruction of the measured nonstationary interference signal by EMD.

Ideally, each IMF contains only one characteristic frequency. However, in the process of signal decomposition, different characteristic frequencies are mixed in the same IMF component, or the same characteristic frequency exists in different IMF components. Wu and Huang etal..^[19] define this phenomenon as modal aliasing. Due to the existence of modal aliasing, the two adjacent IMF components obtained by EMD will interfere with each other, and the characteristic frequency mixing cannot be distinguished. Therefore, in the process of decomposition and reconstruction of non-stationary signals, it is inevitable to eliminate the true component of the signal itself or introduce the false component into the reconstructed signal, which seriously affects the quality of signal reconstruction. As $IM{F_4}$ shown in Figure 3, we generally regard this component as the most important component reflecting the characteristics of the interference signal. Obviously, due to the influence of modal aliasing, other frequency components are mixed into the component $IM{F_4}$, which leads to certain fluctuations in the time domain. Therefore, it directly affects the quality of the reconstructed interference signal ${I_r}(t)$, and then affects the phase extracting accuracy of the interference signal.

In order to overcome the modal aliasing problem, YEH J R etal..^[20] proposed the CEEMD algorithm based on EMD. The CEEMD algorithm successively adds pairs (one positive and one negative) of white noise to the original signal, and then performs EMD separately. The algorithm flow is shown in Figure 4. The CEEMD algorithm uses the characteristics of uniform distribution of white noise in the time-frequency space, and the component signal $IM{F_i}$ of different time-frequency scales is automatically mapped to the appropriate time-frequency scale related to it through the white noise background, to ensure the continuity of each component in the time domain and reduce modal aliasing. In addition, using the time-domain zero mean characteristic of white noise, the CEEMD algorithm eliminates the influence of the added white noise on the original signal by adding white noise pairs multiple times and averaging multiple groups of IMF components of the same order.

The following is the detailed process of CEEMD algorithm decomposition and reconstruction. First, add white noise signal $n(t)$ in pairs to the original signal $s(t)$, that is:

$ \left[ \begin{gathered} {s_{j1}}(t) \\ {s_{j2}}(t) \\ \end{gathered} \right] = \left[ {\begin{array}{*{20}{c}} 1&1 \\ 1&{ - 1} \end{array}} \right]\left[ \begin{gathered} s(t) \\ {n_j}(t) \\ \end{gathered} \right]\quad. $ (22)

The synthesized signals of the original signal $s(t)$ after the $j$th addition of white noise pairs are ${s_{j1}}(t)$ and ${s_{j2}}(t)$, respectively, and EMD of ${s_{j1}}(t)$ and ${s_{j2}}(t)$ will respectively obtain two sets of $n$-layer $IMF$ components, namely:

$ \left[ \begin{gathered} {s_{j1}}(t) \\ {s_{j2}}(t) \\ \end{gathered} \right] = {\left[ {\begin{array}{*{20}{c}} {IM{F_{j1,1}}}&{IM{F_{j2,1}}} \\ {IM{F_{j1,2}}}&{IM{F_{j2,2}}} \\ \vdots & \vdots \\ {IM{F_{j1,n}}}&{IM{F_{j2,n}}} \\ {{r_{j1,n}}(t)}&{{r_{j2,n}}(t)} \end{array}} \right]^{\rm{T}}} \quad.$ (23)

Averaging components of the same order of each layer, we can get:

$ {s_j}(t) = \left[ \begin{array}{*{20}{c}} \dfrac{{IM{F_{j1,1}} + IM{F_{j2,1}}}}{2} \\ \dfrac{{IM{F_{j1,2}} + IM{F_{j2,2}}}}{2} \\ \vdots \\ \dfrac{{IM{F_{j1,n}} + IM{F_{j2,n}}}}{2} \\ \end{array} \right] \quad.$ (24)

When, in Eqs. (22), (23) and (24), $j = 1,2, \cdots , m$, the component matrix of row n and column m can be obtained, that is:

$ \left[ \begin{gathered} {s_1} \\ {s_2} \\ \; \vdots \\ {s_j} \\ \end{gathered} \right] = {\left[ {\begin{array}{*{20}{c}} {IM{F_{1,1}}}&{IM{F_{2,1}}}& \cdots &{IM{F_{m,1}}} \\ {IM{F_{1,2}}}&{IM{F_{2,2}}}& \cdots &{IM{F_{m,2}}} \\ \vdots & \vdots & \ddots & \vdots \\ {IM{F_{1,n}}}&{IM{F_{2,n}}}& \cdots &{IM{F_{m,n}}} \end{array}} \right]^{\rm{T}}}\quad. $ (25)

The component matrix of the CEEMD decomposition signal can be obtained by averaging components of the same order of each layer in the component matrix, namely:

$ s(t) = \left[ \begin{gathered} \frac{1}{m} \cdot \sum\limits_{j = 1}^m {IM{F_{j,1}}} \\ \frac{1}{m} \cdot \sum\limits_{j = 1}^m {IM{F_{j,2}}} \\ \quad\quad\;\;\;\; \vdots \\ \frac{1}{m} \cdot \sum\limits_{j = 1}^m {IM{F_{j,n}}} \\ \end{gathered} \right] = \left[ \begin{gathered} IM{F_1} \\ IM{F_2} \\ \quad \vdots \\ IM{F_n} \\ \end{gathered} \right] \quad.$ (26)

Fig. 5 shows the decomposition and reconstruction results of the measured interference signal by CEEMD algorithm. The components of each order are arranged in order from high frequency to low frequency. By removing the high frequency false components and residual components, the remaining real components are reconstructed to obtain the reconstructed interference signal ${I_{\rm{r}}}(t)$ with high SNR.

Compared with the EMD algorithm, CEEMD algorithm decomposes more components and distinguishes the characteristics between components more clearly, which indicates that CEEMD algorithm can reduce the mutual interference between adjacent IMFs to a certain extent. Comparing the quality of the reconstructed interference signal ${I_{\rm{r}}}(t)$ in Fig. 3 and Fig. 5, the reconstruction quality of the interference signal by the CEEMD algorithm is better than that by the EMD algorithm. Obviously, CEEMD algorithm will help improve the phase solving accuracy of the interference signal. Next, we will verify the effectiveness of the CEEMD algorithm through simulation and experiments.

The FSI ranging platform and experimental related equipment constructed in this experiment are shown in Fig. 6. ECDL (TLB-6712, Newport) and data acquisition card (PXIe-5105, National Instruments) are controlled by the host computer, the frequency modulation range of ECDL is set to 765−775 nm, and the frequency modulation speed is set to 2 nm/s. The optical fiber splitter (TN785R3A1, Thorlabs) divides the light source into two channels, which are coupled to the free space through the collimator (F230APC-780, Thorlabs), and 25% of the light enters the FP etalon (SA200-5B, Thorlabs). The FP signal is detected by the built-in photodetector of the FP etalon, and filtered and amplified by the FP controller; 75% of the light enters the measuring optical path. After passing through the Michelson interference optical path, interference is formed at PD (Model 1801, Newport), and the PD converts the optical interference signal into a voltage signal. DAQ collects FP signal and interference signal synchronously, and the acquired signal will be processed by PC to complete ranging.

To verify the effectiveness of the algorithm, in this paper, the peak point of the transmission signal obtained from the F-P etalon is fitted to obtain the actual optical frequency change curve $\upsilon (t)$ of ECDL output. The optical frequency change rate $\;\beta $ of ECDL output can be expressed as:

$ \beta (t) = \frac{{{\rm{d}}\upsilon (t)}}{{{\rm{d}}t}}\quad. $ (27)

Substituting equation (27) into Eq. (2), and adding noise ${n_p}(t)$ to the simulated interference signal, the final form of the simulated interference signal is obtained as:

$ I'(t) = a \cdot \cos \left[ {2{\text{π}} ({\upsilon _0} + \beta (t) \cdot t) \cdot \frac{{2n \cdot L}}{c}} \right] + b \cdot {n_p}(t)\quad, $ (28)

where a and b are amplitudes of the simulated interference signal and noise respectively.

The simulation parameter settings are shown in Table 1.

In this section, we will compare and verify the ranging performance of CEEMD-HT algorithm and EMD-HT algorithm based on the analysis in Chapter 2 as the theoretical basis. First, we extracted the single cycle phase difference of the interference signal under different signal-to-noise conditions for 20 times and compared the extracted results with the standard single cycle phase difference to reflect the extraction capability of the above two algorithms for decimal cycle phase. As shown in Fig. 7 (color online), both algorithms can reduce the extraction error and standard deviation of the single cycle phase difference compared with the direct extraction of the single cycle phase difference, and under different signal-to-noise conditions, the extract effect of the single cycle phase difference by the CEEMD-HT algorithm is better than that by the EMD-HT algorithm.

Next, to verify the overall measurement performance of the CEEMD-HT algorithm, we have applied the CEEMD-HT algorithm and the EMD-HT algorithm to perform 500 repeated rangings under different SNR conditions and compared them with the standard distance of 10 m set by simulation. The histograms of the measurement residuals and the measurement standard deviation obtained under different SNR conditions are shown in Fig. 8 (color online). The results show that both the CEEMD-HT algorithm and the EMD-HT algorithm can improve the ranging performance of the FSI to some extent compared with the direct measurement.

It is worth noting that when the SNR of the simulated interference signal is 20 dB, the EMD-HT algorithm does not show significant measurement advantages compared with direct measurement. Combined with Fig. 3 and the principle of EMD decomposition reconstruction, it can be seen that EMD depends on the extreme point distribution of the interference signal when decomposing the interference signal, and when the SNR of the simulated interference signal is relatively increased, the volatility of the extreme point distribution decreases, which will aggravate the modal aliasing and have a certain impact on the decomposition performance of the EMD algorithm. Therefore, when the interference signal is reconstructed, some high-frequency noise components will be mixed into the reconstruction signal, and some of the real components of the interference signal will be lost, which will then have a certain impact on the extraction of the phase difference. In addition, combined with the analysis in Chapter 2, when dividing the phase segment, it also depends on the extreme point distribution of the interference signal, so the smoothness of the reconstruction interference signal in the time domain will have a corresponding impact on the interception of the signal segment, which in turn affects the stability of the extraction of the phase difference, resulting in the relative decrease of the measurement stability of the EMD algorithm under this the 20 dB SNR.

Compared with the EMD algorithm, during the decomposition of the interference signal, CEEMD adds white noise pairs to the interference signal to be decomposed, which makes the fluctuation degree of the interference signal with different SNRs in the time domain tend to be consistent. Therefore, the decomposition and reconstruction performance of the CEEMD algorithm is relatively less affected by the SNR. From the simulated ranging results shown in Fig. 8, it can be seen that the measurement residuals and standard deviations of the CEEMD-HT algorithm are lower than those of the EMD-HT algorithm and direct measurements under a certain interference signal SNR. It is worth noting that the CEEMD-HT algorithm can improve the repeatability of the system measurement more effectively than to improve the accuracy of the system measurement. The standard deviation of the CEEMD-HT measurement is 3.14 μm at a SNR of 25 dB for the interference signal and a range of 10 m.

In addition, CEEMD needs to add noise and perform EMD for many times when processing interference signals, which increases the complexity of the algorithm and reduces the efficiency of signal processing to a certain extent. In view of this, the measurement time and efficiency of the CEEMD-HT algorithm are evaluated and analyzed. As shown in Fig. 9 (color online), the measurement efficiency of the CEEMD-HT algorithm is lower than that of the EMD-HT algorithm and direct measurement under the same measurement conditions.

In summary, the CEEMD-HT algorithm can improve the signal processing performance while its simplicity and processing efficiency are also affected accordingly, and it cannot simultaneously take into account the stability and rapidity of the system measurement. However, as an adaptive signal processing means, it can still improve the stability and anti-interference capability of the FSI system without increasing the complexity and cost of the system, which has certain positive significance to improve the performance of the FSI system.

The FSI ranging experimental system is shown in Fig. 6. The sampling rate of the data acquisition device is set to 10 MS/s, and the single sampling time is 5 s, and the FP signal and the interference signal are acquired simultaneously. After the phase segment of the interference signal is divided, the CEEMD-HT algorithm and EMD-HT algorithm are applied to decompose and reconstruct the interference signal of the decimal phase segment respectively, the phase difference of the reconstructed signal is solved, and finally the optical frequency variation range $\Delta \upsilon $ and the phase difference of the intercepted interference signal $\Delta \varphi $ are substituted into the ranging Eq. (7) for ranging.

Before the ranging, we analyzed and compared the performance of CEEMD-HT algorithm and EMD-HT algorithm to process the measured interference signal, and the results are shown in Fig. 10. The measured interference signal disturbed by noise is shown in Fig. 10 (a), and the reconstructed interference signals shown in Fig. 10 (b) and Fig. 10 (c) are obtained after processing by CEEMD algorithm and EMD algorithm. Comparing the locally amplified signals, both the CEEMD algorithm and the EMD algorithm can effectively reduce the interference of noise on the interference signal. However, the interference signal after EMD processing still has fluctuations in the time domain, and such fluctuations will have certain effects on both of performing phase segment division and using HT for phase solution. It can be seen that the CEEMD algorithm is more effective in processing the interference signal under the same conditions.

Based on the above analysis, the actual measurement performance of the CEEMD-HT algorithm is experimentally verified in this section. 50 repetitive measurements are performed on the same target in the 2 m free-space range, and the ranging results are shown in Figure 11. The experimental results show that, in the measurement range of 2 m free space, the repetitive measurement standard deviation based on the CEEMD-HT algorithm is 2.79 μm, which is 5.19 times lower than the repetitive measurement standard deviation of 14.48 μm based on the EMD-HT algorithm, and is 8.28 times lower than the repetitive measurement standard deviation of 23.09 μm by the direct measurement. This shows that the CEEMD-HT algorithm can effectively improve the stability of the FSI system ranging and the anti-interference capability of the system.

Based on the analysis of the effect of optical frequency scanning nonlinearity on the ranging accuracy of the FSI system, we propose a CEEMD-HT-based interference signal phase solving method, and conducts theoretical derivation and simulation analysis on the application feasibility of the CEEMD-HT algorithm in interference signal phase solving. To verify the effectiveness of the CEEMD-HT algorithm for phase solution of non-stationary interference signals, simulation experiments are conducted using the real output optical frequency in the FSI experimental system as the simulation conditions. The simulation results show that the performance of the CEEMD-HT algorithm outperforms that of the direct measurement and EMD-HT algorithm for both single-cycle phase extracting and overall measurement. The standard deviation of the CEEMD-HT algorithm is 3.14 μm in the 10 m range. The measurement efficiency of the algorithms is compared, and the results show that the CEEMD-HT algorithm is inferior to the direct measurement and EMD-HT algorithms in terms of measurement efficiency, and it cannot take into account both the measurement stability and measurement efficiency of the system. In addition, the FSI ranging system was used to verify the measurement results, which showed that, in the measurement range of 2 m free space, the 2.79 μm standard deviation obtained by the CEEMD-HT algorithm is 5.19 times and 8.28 times lower than that of the EMD-HT algorithm and the direct measurement method, respectively. Therefore, the CEEMD-HT algorithm can effectively reduce the effects of optical frequency scanning nonlinearity and noise on the phase extracting and ranging stability of the interference signal, and improve the anti-interference and robustness of the FSI ranging system.

光频扫描干涉（Frequency Scanning Interferometry，FSI）绝对测距技术在多卫星编队、大型装备装配制造、自动驾驶等领域有着广泛的应用前景^[1-9]。FSI测距系统的典型光源为窄线宽可调谐外腔半导体激光器(External Cavity Diode Laser，ECDL)，其通过压电陶瓷或电机改变外腔长度实现对光频率的连续调谐。然而，受到位移调谐元件非线性响应的影响，在线性驱动信号输入下，ECDL输出的光频率变化与时间呈非线性关系^[10-12]。ECDL的光频扫描非线性导致FSI干涉信号的频率时变特性，进而影响了干涉信号相位的提取精度，最终导致了FSI测距精度下降^{[2, 11, 13-14]}。因此，降低光频扫描非线性对干涉信号的影响对于提高FSI测距精度至关重要。

当前针对ECDL的光频扫描非线性问题的研究方法可归纳为以下两类：第一类方法是通过监测光频变化信息，并利用反馈环路对光频扫描速率进行在线控制，实现扫频非线性的前期抑制，如2009年，Roos P A等人^[15]利用双锁相环法，对ECDL的光频描速率进行在线控制, 降低了扫频非线性对干涉信号频率的影响，在1.5 m范围内实现了86 nm的测距精度。这类方法通过引入反馈控制环路，在光路形成干涉前对扫频非线性进行抑制，然而引入的反馈控制环路增加了系统复杂度，环路带宽同时也限制了光频扫描速率的提高；另一类方法则通过对频率时变的干涉信号采样后进行后处理，消除扫频非线性对干涉信号相位提取精度的影响，如2014年，时光等人^[2]构造了30 m光程差的全光纤马赫-曾德参考干涉仪，利用生成的参考干涉信号对测量干涉信号进行等光频间隔重采样，大幅降低了扫频非线性对测量干涉信号的影响，将测距系统空间分辨率提高了约30倍。然而在此类重采样方法中，用于重采样的参考干涉信号需要满足奈奎斯特采样定理，导致系统测距范围受到参考干涉仪光程差的限制，此外参考干涉仪的引入同样增加了测距系统的复杂度。

以上两类方法都不同程度地增加了系统复杂度，并且限制了测距系统的测量带宽和测量范围等性能。直到2018年，邓文等人^[16-18]提出，使用经验模式分解(Empirical Mode Decomposition，EMD)结合希尔伯特变换（Hilber Transform，HT）对非平稳干涉信号进行相位提取，该方法有效地降低了光频扫描非线性对干涉信号相位提取精度的影响。然而，传统EMD在对非平稳信号进行分解重构时存在分量模态混叠现象，限制了相位提取精度的进一步提高。为此，本文利用基于噪声辅助分析的非平稳信号模式分解重构方法—互补集合经验模态分解( Complementary Ensemble Empirical Mode Decomposition，CEEMD)对测量干涉信号进行分解重构，改善传统EMD方法造成的模态混叠现象，进一步降低光频扫描非线性的影响。

如图1所示，FSI测距系统由ECDL输出扫频激光，经准直后在分光镜(BS1)处被分成两束：一束进入法布里-珀罗（Fabry-Perot, F-P）标准具，经光电探测器PD1探测产生透射F-P信号；另一束进入迈克尔逊干涉光路，并由分光镜(BS2)分为参考光和测量光，分别经反射棱镜RR1和RR2反射后在光电探测器PD2上形成干涉信号。

在FSI测距系统中，ECDL的输出光频率 $\upsilon (t) $可表示为：

$ \upsilon (t)={\upsilon }_{0}+\beta \cdot t\quad, $ (1)

式中： $ {\upsilon }_{0} $为激光器输出的初始光频率， $ \;\beta $为光频扫描速率。根据光波干涉原理，干涉信号可表示为：

$ I(t)={I}_{{\rm{r}}}+{I}_{{\rm{m}}}+2\sqrt{{I}_{{\rm{r}}}\cdot {I}_{{\rm{m}}}}\mathrm{cos}\left[2{\text{π}} {\upsilon }_{0}\tau +2{\text{π}} \beta \cdot \tau \cdot t\right]\quad, $ (2)

式中， $ {I}_{{\rm{r}}} $为参考光光强， $ {I}_{{\rm{m}}} $为测量光光强， $ \tau $为测量光相对于参考光在时域上的延迟。由式(2)可得，干涉信号的瞬时相位 $ \varphi $表示为：

$ \varphi = 2{\text{π}} {\upsilon _0}\tau + 2{\text{π}} \beta \cdot \tau \cdot t \quad. $ (3)

若在 $ \Delta t $=t₂−t₁时间内，ECDL的光频变化量为 $ \Delta \upsilon $，则对应的相位差 $ \Delta \varphi $可表示为：

$ \Delta \varphi = \varphi ({t_2}) - \varphi ({t_1}) = 2{\text{π}} \cdot \tau \cdot \beta \cdot \Delta t \quad,$ (4)

结合式(1)，相位差 $ \Delta \varphi $可进一步改写为：

$ \Delta \varphi = 2{\text{π}} \cdot \tau \cdot \left[ {\upsilon ({t_2}) - \upsilon ({t_1})} \right] = 2{\text{π}} \cdot \tau \cdot \Delta \upsilon\quad.$ (5)

当待测距离为 $ L $时，延时 $ \tau $可表示为：

$ \tau = \frac{{2n \cdot L}}{c} \quad,$ (6)

式中， $ n $为空气折射率， $ c $为光速。联立式(5)、(6)可得待测距离L的表达式：

$ L = \frac{{c \cdot \tau }}{{2n}} = \frac{{c \cdot \Delta \varphi }}{{4{\text{π}} \cdot n \cdot \Delta \upsilon }}\quad. $ (7)

由式(7)可知，系统的测距精度取决于干涉信号相位差 $ \Delta \varphi $的提取和光频变化范围 $ \Delta \upsilon $的计算。光频变化范围可通过对图1中的F-P信号波峰数量计数获得，选定起始、终止F-P波峰后，光频变化范围 $ \Delta \upsilon $与F-P波峰数 $ r $间存在如下关系：

$ \Delta \upsilon = (r - 1) \cdot {FSR}\quad, $ (8)

式中， $ \text{FSR} $为F-P标准具的自由光谱范围（Free Spectrum Range, FSR），即相邻峰之间对应的光频差。因此，求解扫频范围 $ \Delta \upsilon $对应的干涉信号相位差 $ \Delta \varphi $是决定FSI测距精度的关键。

事实上，经FP信号截取后的干涉信号包含了整数周期信号段和小数周期信号段，整数周期信号段内相邻两个波峰间对应的相位差为 $ 2{\text{π}} $，因此整数周期相位差可直接通过干涉信号的波峰数求解。如图2所示，为提高干涉信号相位差的求解效率，将干涉信号的相位差分为三部分求解，整数相位差 $ \Delta {\varphi }_{{\rm{i}}} $、起始小数相位差 $ \Delta {\varphi }_{{\rm{s}}} $和截止小数相位差 $ \Delta {\varphi }_{{\rm{e}}} $。在进行相位段自动截取时，以起始F-P波峰所在时刻 $ {t}_{{\rm{s}}} $为基准向右索引干涉的第一个波峰所在时刻 $ {t}_{{\rm{p}}1} $，该区间对应为起始小数相位差区间；同理，以终止F-P波峰所在时刻 $ {t}_{{\rm{e}}} $为基准向左索引干涉的最后一个波峰所在时刻 $ {t}_{{\rm{p}}2} $，该区间则为截止小数相位区间；而 $ {t}_{{\rm{p}}1} $到 $ {t}_{{\rm{p}}2} $所对应的区间即为整数相位区间。整数相位差 $ \Delta {\varphi }_{{\rm{i}}} $可表示为：

$ \Delta {\varphi _{\rm{i}}} = (\delta - 1) \cdot 2{\text{π}}\quad. $ (9)

整数相位差由干涉信号的峰值数 $ \delta $决定，小数相位差 $ \Delta {\varphi }_{{\rm{s}}} $， $ \Delta {\varphi }_{{\rm{e}}} $则利用HT进行求解：

$ \left\{ \begin{gathered} \Delta {\varphi _{\rm{s}}} = {\varphi _{\rm{u}}}({t_{{\rm{p}}1}} - {t_{\rm{s}}})/T \\ \Delta {\varphi _{\rm{e}}} = {\varphi _{\rm{u}}}({t_{\rm{e}}} - {t_{{\rm{p}}2}})/T \\ \end{gathered} \right. \quad, $ (10)

上式中， $ {\varphi }_{{\rm{u}}} $为经HT后的相位解包裹函数， $ T $为干涉信号的周期。因此，干涉信号的相位可表示为：

$ \Delta \varphi = \Delta {\varphi _{\rm{i}}} + \Delta {\varphi _{\rm{s}}} + \Delta {\varphi _{\rm{e}}} \quad. $ (11)

此时，将扫频范围 $ \Delta \upsilon $和对应的干涉信号相位差 $ \Delta \varphi $代入式(7)即可解算得到被测距离 $ L $。

然而，ECDL的位移调谐元件PZT存在迟滞和蠕变特性，导致在线性驱动信号输入下，ECDL输出的光频率变化与时间呈非线性关系，光频扫描速率 $ \;\beta $不再为定值，干涉信号周期T成为时变量，利用HT计算小数相位时将会产生误差。因此，本文提出利用CEEMD-HT算法计算所截干涉信号段内的逐点小数相位值，进而实现所截干涉信号相位变化量的计算，以避免扫频非线性对其相位计算的影响。

首先对HT在光频扫描干涉信号相位求解方法进行推导，对于任意信号 $ x(t) $，其希尔伯特变换可表示为：

$ H[x(t)]={\mathcal{F}}^-[\mathcal{F}[x(t)]\delta (\omega )] \quad,$ (12)

式(12)中 $\mathcal{F}[\cdot ]$和 ${\mathcal{F}}^{-}[\cdot ]$为傅立叶变换和傅立叶逆变换， $ \delta (\omega ) $为HT的窗函数，其可表示为：

$\delta (\omega ) = - j{{\rm{sgn}}} (\omega ) = \left\{ \begin{gathered} - j{\text{ }}(w > 0) \\ 0{\text{ (}}w = 0{\text{)}} \\ j{\text{ (}}w < 0{\text{)}} \\ \end{gathered} \right.\quad. $ (13)

由式(2)可知，FSI干涉信号 $ I(t) $表达式可简化为下式：

$ I(t) = a(t) + b(t)\cos \varphi (t)\quad, $ (14)

式中， $ a(t) $为干涉信号的直流分量， $ b(t) $为干涉信号的幅度调制，干涉信号的相位信息隐含于 $ b(t)\mathrm{cos}\varphi (t) $中。利用合适的高通器滤波掉低频直流项 $ a(t) $，得到包含干涉信号相位信息的 $ b(t)\mathrm{cos}\varphi (t) $项，此时干涉信号可表示为下式：

$ {I_p}(t) = {\mathcal{F}^ - }[\mathcal{F}[I(t)] \cdot W(\omega )] = b(t)\cos \varphi (t)\quad,$ (15)

式(15)中， $ W(\omega ) $为高通滤波器的窗函数。此外，Bedrosian E^[17]等人对乘积形式的希尔伯特变换做了详细的论述，并得到了重要结论：低通信号与高通信号乘积的希尔伯特变换等价于低通信号与高通信号希尔伯特变换的乘积。基于此结论，高通滤波后的干涉信号的希尔伯特变换可表示为：

$ H[{I_p}(t)] = b(t)H[\cos \varphi (t)] = b(t)\sin \varphi (t) \quad. $ (16)

在上述分析的基础上，利用希尔伯特变换构造干涉信号的解析函数 $ \hat{I}(t) $：

$ \hat I(t) = {I_p}(t) + iH[{I_p}(t)]\quad. $ (17)

结合式(15)，式(16)，可得：

$ \hat I(t) = b(t)\cos \varphi (t) + ib(t)\sin \varphi (t)\quad. $ (18)

此时，干涉信号的相位函数 $ \varphi (t) $可表示为：

$ \varphi (t) = \arctan \left[ {\frac{{{{\rm{Im}}} (\hat I(t))}}{{{{\rm{Re}}} (\hat I(t))}}} \right] = \arctan \left[ {\frac{{b(t)\sin \varphi (t)}}{{b(t)\cos \varphi (t)}}} \right]\quad. $ (19)

由式(19)可得， $ {t}_{1} $时刻到 $ {t}_{2} $时刻，干涉信号 $ I(t) $的相位变化量 $ \Delta \varphi $可表示为：

$ \Delta \varphi = {\varphi _{\rm{u}}}({t_2}) - {\varphi _{\rm{u}}}({t_1})\quad. $ (20)

需要注意的是，干涉信号的信噪比是影响HT相位求解精度的关键因素^[16]。干涉信号受测量环境扰动、光电探测器噪声等因素影响，尤其在长距离测量条件下，其信噪比往往较低，故此处利用EMD对干涉信号进行自适应分解，并重构出高质量的干涉信号，提高干涉信号相位提取精度。在信号分解方面EMD并不依赖于基函数，完全根据信号本身的极值点分布进行自适应分解，将包含多种频率成分的复杂信号逐级分解为多个单频分量之和的形式，适用于具有非平稳特性的复杂信号的分析处理，EMD分解的具体过程描述如下:

(1) 先提取原始信号 $ s(t) $的上下极值点, 分别用三次样条拟合上下极值点得到上包络信号 $ {e}_{{\rm{u}}}(t) $和下包络信号 $ {e}_{{\rm{d}}}(t) $。

(2) 求得上下包络的均值包络信号 $ m(t) $，从初始信号 $ s(t) $中剔除掉 $ m(t) $，得到信号 $ c(t) $。

(3) 将 $ c(t) $视为新的 $ s(t) $信号，重复步骤(1)~(2)，直到最新分解出的 $ c(t) $满足分解终止条件^[17]，此时 $ c(t) $视为第一阶分量 $ IM{F}_{1} $。

(4) 从 $ s(t) $中剔除掉 $ IM{F}_{1} $，剩余信号记作 $ r(t) $，将 $ r(t) $视为新的 $ s(t) $信号，重复步骤(1)~(3)，依次得到各阶分量 $ IM{F}_{i} $。

非平稳信号 $ s(t) $的分解形式可表示为：

$ s(t) = \sum\limits_{i = 1}^n {IM{F_i}} + {r_n}(t)\quad,$ (21)

式中， $ IM{F}_{1} $， $ IM{F}_{2},\cdots, $$ IM{F}_{n} $为EMD分解所得的各阶分量，也称为本征模态函数(Intrinsic Modal Function，IMF)，图3为EMD对实测非平稳干涉信号的分解重构示意图。

在理想情况下，每个IMF中仅包含一种特征频率。然而，EMD在对信号分解的过程中，同一个IMF分量中混杂着不同的特征频率，或是不同的IMF中存在着相同的特征频率，WU ZH H和 HUANG N E等^[19]将这种现象定义为模态混叠。由于模态混叠的存在，导致EMD分解得到的相邻两个IMF分量之间存在相互干扰，特征频率混杂，无法区分。因此，在对非平稳信号进行分解重构的过程中，不可避免地会剔除掉信号本身的真实分量或是将虚假分量引入到重构信号中，以至于严重影响了信号的重构质量。如图3中的 $ IM{F}_{4} $，一般将该分量视为反映干涉信号特征的最主要的分量。显然，受模态混叠的影响， $ IM{F}_{4} $分量中混入了其他频率成分，导致其在时域上产生了一定的波动，因此直接影响了重构干涉信号 $ {I}_{r}(t) $的质量，继而影响了干涉信号的相位提取精度。

为了克服模态混叠问题，YEH J R等^[20]在EMD的基础上提出了CEEMD算法。CEEMD算法先后向原始信号中添加成对(一正一负)的白噪声，再分别进行EMD分解，算法流程如图4所示。CEEMD算法利用白噪声在时频空间上均匀分布的特性，将不同时频尺度的分量信号 $ IM{F}_{i} $通过白噪声背景自动映射到与之相关的适当的时频尺度上，以保证每个分量在时域上的连续性进而减少模态混叠。此外，利用白噪声的时域零均值特性，CEEMD算法通过多次添加白噪声对，并对多组同阶的IMF分量取平均，以消除所加白噪声对原始信号的影响。

以下是CEEMD算法分解重构的详细过程，先向原始信号 $ s(t) $中添加成对的白噪声信号 $ n(t) $，即：

$ \left[ \begin{gathered} {s_{j1}}(t) \\ {s_{j2}}(t) \\ \end{gathered} \right] = \left[ {\begin{array}{*{20}{c}} 1&1 \\ 1&{ - 1} \end{array}} \right]\left[ \begin{gathered} s(t) \\ {n_j}(t) \\ \end{gathered} \right]\quad, $ (22)

原始信号 $ s(t) $第 $ j $次添加白噪声对后的合成信号分别为 $ {s}_{j1}(t) $和 $ {s}_{j2}(t) $，分别对 $ {s}_{j1}(t) $和 $ {s}_{j2}(t) $进行EMD分解，每次分解将得到两组 $ n $层 $ IMF $分量，即：

$ \left[ \begin{gathered} {s_{j1}}(t) \\ {s_{j2}}(t) \\ \end{gathered} \right] = {\left[ {\begin{array}{*{20}{c}} {IM{F_{j1,1}}}&{IM{F_{j2,1}}} \\ {IM{F_{j1,2}}}&{IM{F_{j2,2}}} \\ \vdots & \vdots \\ {IM{F_{j1,n}}}&{IM{F_{j2,n}}} \\ {{r_{j1,n}}(t)}&{{r_{j2,n}}(t)} \end{array}} \right]^{\rm{T}}} \quad. $ (23)

对各层的同阶分量取平均，可得：

$ {{\boldsymbol{s}}_j}(t) = \left[ \begin{array}{*{20}{c}} \dfrac{{IM{F_{j1,1}} + IM{F_{j2,1}}}}{1} \\ \dfrac{{IM{F_{j1,2}} + IM{F_{j2,2}}}}{2} \\ \vdots \\ \dfrac{{IM{F_{j1,n}} + IM{F_{j2,n}}}}{2} \\ \end{array} \right] \quad. $ (24)

当在式(22)~(24)中， $ j=1,2,\cdots ,m $时，可得到 $ n $行 $ m $列的分量矩阵，即：

$ \left[ \begin{gathered} {s_1} \\ {s_2} \\ \; \vdots \\ {s_j} \\ \end{gathered} \right] = {\left[ {\begin{array}{*{20}{c}} {IM{F_{1,1}}}&{IM{F_{2,1}}}& \cdots &{IM{F_{j,1}}} \\ {IM{F_{1,2}}}&{IM{F_{2,2}}}& \cdots &{IM{F_{j,2}}} \\ \vdots & \vdots & \ddots & \vdots \\ {IM{F_{1,n}}}&{IM{F_{2,n}}}& \cdots &{IM{F_{j,n}}} \end{array}} \right]^{\rm{T}}}\quad. $ (25)

对分量矩阵中的各层同阶分量求平均，可得CEEMD分解信号的分量矩阵，即：

$ {\boldsymbol{s}}(t) = \left[ \begin{gathered} \frac{1}{m} \cdot \sum\limits_{j = 1}^m {IM{F_{j,1}}} \\ \frac{1}{m} \cdot \sum\limits_{j = 1}^m {IM{F_{j,2}}} \\ \quad\quad\;\;\;\; \vdots \\ \frac{1}{m} \cdot \sum\limits_{j = 1}^m {IM{F_{j,n}}} \\ \end{gathered} \right] = \left[ \begin{gathered} IM{F_1} \\ IM{F_2} \\ \quad \vdots \\ IM{F_n} \\ \end{gathered} \right] \quad. $ (26)

图5为CEEMD算法对实测干涉信号的分解重构示意图，各阶分量按频率从大到小依次排列，通过去除掉高频虚假分量和残余分量，对剩余的真实分量进行重构，得到高信噪比的重构干涉信号 $ {I}_{{\rm{r}}}(t) $。

与EMD算法分解相比，CEEMD算法分解出的分量更多，各分量间的特征区分更为清晰，这表明CEEMD算法在一定程度上能够降低相邻IMF之间的相互干扰。对比图3和图5中重构干涉信号 $ {I}_{{\rm{r}}}(t) $的质量可知，CEEMD算法对干涉信号的重构质量要优于EMD算法，显然，CEEMD算法将有助于提高干涉信号的相位求解精度。接下来，将通过仿真和实验对CEEMD算法的有效性进行验证。

本实验构建的FSI测距平台及实验相关设备如图6所示，通过上位机对ECDL(TLB-6712, Newport)和数据采集卡(PXIe-5105, National Instruments)进行控制，设定ECDL的调频范围为765~775 nm，调频速度为2 nm/s。光纤分路器(TN785R3A1, Thorlabs)将光源分为两路，经准直器(F230APC-780，Thorlabs)耦合到自由空间，其中，25%的光进入FP标准具(SA200-5B, Thorlabs)，由FP标准具内置光电探测器检测FP信号，经FP控制器将FP信号滤波放大；75%的光进入测量光路，经过迈克尔逊干涉光路后，在PD(Model 1801, Newport)处形成干涉，由PD将光干涉信号转化为电压信号。DAQ对FP信号和干涉信号进行同步采集，PC端将对采集到的信号进行处理，完成距离测量。

为验证算法的有效性，本文对由F-P标准具所获得的透射信号的峰值点进行拟合，得到ECDL输出的实际光频变化曲线 $ \upsilon (t) $，此时ECDL输出的光频扫描速率 $ \;\beta $可表示为：

$ \beta (t) = \frac{{{\rm{d}}\upsilon (t)}}{{{\rm{d}}t}}\quad, $ (27)

将式(27)代入式(2)，并向仿真干涉信号中添加噪声 $ {n}_{p}(t) $，得到仿真干涉信号的最终形式：

$ I'(t) = a \cdot \cos \left[ {2{\text{π}} ({\upsilon _0} + \beta (t) \cdot t) \cdot \frac{{2n \cdot L}}{c}} \right] + b \cdot {n_p}(t)\quad, $ (28)

式中，a为干涉信号的幅值，b为噪声信号幅值。仿真参数设置如表1所示。

本节，将以第2章的分析为理论依据，对CEEMD-HT算法和EMD-HT算法的测距性能进行对比验证。首先，对不同信噪比下的干涉信号的单周期相位差进行了20次重复提取，并将提取结果与标准的单周期相位差进行比对，以此来反映上述两种算法对小数周期相位的提取能力。如图7所示，与直接提取的单周期相位差相比，两种算法均能够降低单周期相位差的提取误差和标准差，同时，不同信噪比条件下，CEEMD-HT算法对单周期相位差的提取效果均优于EMD-HT算法。

接下来，为验证CEEMD-HT算法的总体测量性能，先后运用CEEMD-HT算法和EMD-HT算法对不同信噪比条件下的距离进行了500次重复测量，并与仿真设置的标准距离10 m进行了对比。如图8所示为不同信噪比条件下得到的测量残差与测量标准差的直方图。结果表明，与直接测量相比，CEEMD-HT算法和EMD-HT算法在一定程度上均能够提高FSI的测距性能。

值得注意的是，在仿真干涉信号信噪比为20 dB时，与直接测量相比，EMD-HT算法并未表现出明显的测量优势。结合图3与EMD分解重构原理进一步分析可知，EMD在对干涉信号进行分解时，需要依赖干涉信号的极值点分布，当仿真干涉信号信噪比相对提高时，其极值点分布的波动性降低，该情况会加剧模态混叠，对EMD算法的分解性能产生一定的影响。因而，在进行干涉信号重构时，不但会使部分高频噪声分量混入到重构信号当中，同时也会造成干涉信号真实分量的损失，继而对相位差的提取产生一定的影响，此外，结合第2章的分析，在进行相位段划分时，也依赖于干涉信号的极值点分布，因此，重构干涉信号在时域上的平滑程度会对信号段的截取产生相应的影响，继而影响到相位差提取的稳定性，导致该信噪比条件下EMD算法测量稳定性的相对下降。

与EMD算法相比，CEEMD在对干涉信号进行分解的过程中，向待分解的干涉信号中加入白噪声对，这使得不同信噪比的干涉信号在时域上的波动程度趋于一致，因此，CEEMD算法的分解重构性能受信噪比的影响相对较小。由图8所示的仿真测距结果可知，在干涉信号信噪比一定的条件下，CEEMD-HT算法的测量残差和标准差均要低于EMD-HT算法和直接测量的结果。值得注意的是，相较于改善系统测量的准确性，CEEMD-HT算法能够更为有效地提高系统测量的重复性。在干涉信号信噪比为25 dB时，10 m测距范围内，CEEMD-HT的测量标准差为3.14 μm。

此外，CEEMD在对干涉信号进行处理时，需要多次添加噪声及进行多次EMD分解，这在一定程度上增加了算法的复杂性，降低了信号处理的效率，而对于测量而言，测量快速性是不容忽视的性能指标。鉴于此，本文对CEEMD-HT算法的测量时间和效率进行了评估分析。如图9所示，在同一测量条件下，CEEMD-HT算法的测量效率要低于EMD-HT算法和直接测量。

综上所述， CEEMD-HT算法在提升信号处理性能的同时，其算法简洁性和处理效率也受到了相应的影响，无法同时兼顾系统测量的稳定性和快速性。但是，作为一种自适应信号处理手段，在不增加系统复杂度和成本的基础上，依旧能够提升FSI系统的稳定性和抗干扰能力，对改善FSI系统的性能具有一定的积极意义。

FSI测距实验系统如图6所示，设置数据采集装置的采样率为10 MS/s，单次采样时间为5 s，并对 FP信号和干涉信号进行同步采集。在对干涉信号相位段进行划分后，分别运用CEEMD-HT算法和EMD-HT算法对小数相位段的干涉信号进行分解重构，并对重构信号的相位差进行求解，最后将光频变化范围 $ \Delta \upsilon $和所截取的干涉信号的相位差 $ \Delta \varphi $代入测距公式(7)，进行距离测量。

在进行距离测量之前，对 CEEMD-HT算法和EMD-HT算法处理实测干涉信号的性能进行了对比分析，结果如图10所示，图10(a)所示为受噪声干扰的实测干涉信号，经CEEMD算法和EMD算法处理后，得到如图10(b)和图10(c)所示的重构干涉信号。对比局部放大信号可知，CEEMD算法和EMD算法均能有效降低噪声对干涉信号的干扰。但是，经EMD处理后的干涉信号在时域上仍然存在波动，这种波动无论是在进行相位段划分时，还是使用HT进行相位求解时均会产生一定的影响。由此可见，相同条件下CEEMD算法对干涉信号的处理效果更为明显。

在上述分析的基础上，本节对CEEMD-HT算法的实际测量性能进行了实验验证，在2 m自由空间范围内对同一目标进行了50次重复测量，测距结果如图11所示。实验结果表明：在2 m自由空间测量范围内，基于CEEMD-HT算法的重复测量标准差为2.79 μm，较EMD-HT算法的重复测量标准差14.48 μm降低了5.19倍，较直接测量重复测量标准差23.09 μm降低了8.28倍。由此可见，CEEMD-HT算法能够有效提高FSI系统测距的稳定性和系统的抗干扰能力。

本文在分析光频扫描非线性对FSI系统测距精度影响的基础上，提出了一种基于CEEMD-HT的干涉信号相位求解方法，并对CEEMD-HT算法在干涉信号相位求解中的应用可行性进行了理论推导和仿真分析。为验证CEEMD-HT算法对非平稳干涉信号相位求解的有效性，采用FSI实验系统中的真实输出光频率作为仿真条件，进行了仿真实验。仿真结果表明：无论是单周期相位提取还是整体测量，CEEMD-HT算法的性能表现均优于直接测量和EMD-HT算法的性能；在10 m仿真测距范围内，CEEMD-HT算法的测距标准差为3.14 μm；同时，对算法的测量效率进行了分析，结果表明：CEEMD-HT算法在测量效率方面要逊色于直接测量和EMD-HT算法，其无法同时兼顾系统的测量稳定性和测量效率。此外，利用FSI测距系统进行了实测验证，测量结果表明：在2 m自由空间测量范围内，基于CEEMD-HT算法的测距标准差为2.79 μm，相较于EMD-HT算法和直接测量法测量标准差分别降低了5.19倍和8.28倍。因此，CEEMD-HT算法能够有效降低光频扫描非线性以及噪声对干涉信号相位提取及测距稳定性的影响，提高FSI测距系统的抗干扰性和鲁棒性。