Multiwavelength Frequency Modulated CW Ladar: The Effect of Refractive Index

: Frequency modulated continuous wave (FMCW) laser detection and ranging is a technique for absolute distance measurements with high performances in terms of resolution, non-ambiguity range, accuracy and fast detection. It is based on a simple experimental setup, thus resulting in cost restraint with potential wide spread, not only limited to research institutions. The technique has been widely studied and improved both in terms of experimental setup by absolute reference or active stabilization and in terms of data analysis. Very recently a multi-wavelength approach has been exploited, demonstrating high precision and non ambiguity range. The variability of refractive index along the path was not taken into account with consequent degradation of range accuracy. In this work we developed a simple model able to take into account refractive index effect in multi-wavelength FMCW measurement. We performed a numerical simulation in different atmospheric conditions of temperature, pressure, humidity and CO 2 concentration showing a net improvement of range accuracy when refractive index modeling is used.


Introduction
Precision measurement of non-cooperative target by laser ranging (LR) techniques had and continues to have applications in several research fields: from satellites flying in defined formation for large synthetic aperture telescope [1][2][3][4] used for fundamental physics test or extraterrestrial planets search to earth based laser ranging. The earth based satellite laser ranging was successfully demonstrated in 1964 and, since then, produced a huge amount of results in fundamental physics experiments like: general relativity test (gravitomagnetic effect [5], Lense-Thirring [6,7] non-Newtonian gravity [8]) or dark matter search [9,10]. However, it mostly finds application in fields like geodesy and geodynamics, allowing , for example, the length-of-day determination or earth gravitational field measurement [11] and contributing to the definition and the update of the terrestrial reference frame (TRF). In addition, there are perspectives in space debris applications and multipurpose communication experiments. The precise distance determination of near object has also several applications mostly in industrial fields [12], particularly in automotive [13,14]. The main characteristics of laser ranging measurements are: resolution of the detected length, non-ambiguity range and accuracy. A LR system is based on a cw or pulsed laser source, collimation optics, a detector and a timing system. A system based on pulsed source has large non-ambiguity range but low resolution whereas a cw system shows a very small non-ambiguity range (of the order of the wavelength) but extraordinary resolution. The accuracy, instead, does not depend on the particular laser source used. First of all, it is related to experimental setup calibration (e.g., linearity of frequency sweep, electrical cable length calibration or accuracy of timing system). Moreover, different environmental factors may influence the accuracy: one of the main contributor is the refractive index. A great advance was the introduction of multi-wavelength interferometry [15,16] that combines measurements at different wavelengths, generating a longer 'synthetic wavelength' and, therefore, a larger ambiguity range while maintaining same resolution. Unfortunately, extending the ambiguity range beyond a millimetre is very hard. In 2009 Coddington et al. [17] used dual coherent frequency combs approach for ranging experiments, obtaining at same time large non-ambiguity range, exceptional resolution and fast refresh rate. On the downside, the setup used is very complex and expensive and so available only in state-of-the-art laboratories. For this reason a very popular method that keeps costs low but good performance is the frequency modulated continuous wave ranging (FMCW) [18][19][20]. In this technique a cw laser is chirped at a constant rate K (Hz/s) and split in two parts: reference and probe beam. The probe beam is collimated and sent to object whose distance must be measured, and the reflection is mixed with reference beam on a photodetector. The electric signal is a beat note at frequency f = 2Lk/c where L is the optical path, and c the speed of light. The performance of FMCW depends on the coherence length of the used laser (related to spectral full width half maximum) and on total laser frequency excursion [21] that, in turn, depends on chirp rate K. The accuracy is related to linearity of chirp rate and to the equation used to model refractive index. The non-linearity of the chirp rate can be compensated in an active way using a real time detection system and actuators [22,23]. Another way is a post-processing correction by using a fine calibration up to 15 ppb accuracy level by using a frequency comb [24][25][26] or ppm level with a simple and cheaper reference against a molecular spectrum [27].
Recently multi-wavelength approach has been applied to FMCW [28,29] demonstrating high precision and non ambiguity range. In these experimental implementations the variability of refractive index is not considered. Consequently, even if the range uncertainty is low, the overall accuracy of range measurement is large since it does not consider the refractive index that masks the absolute range value introducing a systematic uncertainty. In this paper we develop a model to take into account the air refractive index in a FMCW multi-wavelength approach. The multiple observables introduced in multi-wavelength approach allow to measure not only the range but also the average atmospheric parameters (like pressure, temperature, humidity and CO 2 concentration) along the path. Such parameters allow to calculate a sort of effective refractive index along the path enabling a more accurate estimation of real absolute distance. In other words, in single wavelength FMCW the atmospheric parameters are measured at laser station and used to calculate the refractive index. In this way the refractive index along the path is approximated with the one at laser station. The model developed in this paper allow to measure the effective refractive index along the path in multiwavelength FMCW. Consequently, the refractive index along the path is approximated with the effective refractive index (and not with index at laser station) with consequent improved accuracy. Finally a comparison between two possible atmospheric conditions and a comparison between possible multiple wavelengths has been performed.

Model
A general equation for object range (R) detection via beat frequency measurement of a chirped signal (linearly-chirped with constant K) is j means that we are applying the same equation for R but using different values of laser wavelength λ j hence measuring a different possibly specific beat frequency f j .
Following the formula proposed by Ciddor [30] for the excess component in air n(x, λ) compared to the refractive index of the vacuum, the former depends on λ j (from 0.3 to 1.69 µm in [30], formula validity extended in mid infrared [31]), t the temperature in Celsius (−40 to +100 • C), p the pressure in Pascal (80 to 120 kPa), h the fractional humidity in the [0, 1] interval, x c the CO 2 concentration (from 0 to 2000 ppm). Please see https://emtoolbox.nist.gov/Wavelength/Documentation.asp for detailed explanation of the Ciddor formula used.

First order Taylor Approximation
Using first order Taylor approximation around a known refractive index value n(x (0) , λ j ), and using the notation of sum over repeated indexes: In case of the Ciddor formula l = 4, as we can see in Equation (2). Please notice that if we know c ij , i.e. x (0) and the derivatives of n, and if we have l + 1 = m, we can invert the equations isolating known and unknown quantities: In the last equation we have inverted the linear system generated by the Taylor expansion of n (C −1 is in fact the inverse of matrix C). In Equation (11) we have redefined the index convention for ease of notation in the following treatment. The elements of C −1 are redefined as a ij , where now j sums up with the beat frequencies in F while i identifies the components of the parameters vector P where the range R and the atmospheric parameters x are stored. Such convention will be kept from now on unless otherwise stated. From the above we can find R and x values in terms of only known parameters: So the variance for R and x i will be (covariances excluded):

Effect of Refractive Index Variability
In Equation (1) we assumed the refractive index has a constant value along the laser trajectory, depending only on the constant values of the atmospheric parameters x. In this section we start relaxing this assumption and provide an estimation of the error due to the refractive index variability.
We can think of the refractive index variation along the light path as an accumulation of d f j in Equation (1) along the way. Namely the longer the path, the bigger will be the beat frequency observed. To test the performance of the proposed approach, we evaluate the integral (16) using simplified models for the variation of atmospheric parameters x between the target and the laboratory where the measure is performed.
In the next section we consider the case of x fluctuating around an average value; this simplifies the analytical treatment and the exact solution, allowing for direct estimation over the entire validity range of Ciddor formula. Then we generalize to the case of a linear gradient of x between lab and target, solving numerically integral (16) also for very extreme differences between lab and target in order to test the performance of the proposed approach in limit cases of input variability. We also sketch up the analytical treatment in the general case when x has an unknown three-dimensional spatial distribution to be chosen.

Refractive Index Fluctuation around an Average Value
In general the atmospheric variables x have a three-dimensional distribution in space. For instance, increasing altitude of the target means the light trajectory passes through a gradient of temperature; same usually happens with pressure and the other parameters x. In this section we will assume that the variability of n is dependent on fluctuations of the x spatial distribution around an average value with a certain distribution g. This allows an easier analytical treatment of (16) producing a simple formula to test the performance of the proposed approach that can be evaluated over the entire x validity range of Ciddor formula. Such an assumption can be mathematically stated as follows: Hence from (16) This last equation gives an estimation of the beat frequency when R is known. We can substitute it in Equation (12)-dropping the first index for simplicity of notation so that now a 1j is a j σ rel R * can be seen as an estimate of the relative error due to the usage of Equation (12) for range estimation R * . Notice that even when E[n j ] goes to zero, such error does not necessary go to zero as well. E[n j ] could be estimated as n(E[x], λ j ); the a j are estimated from the derivatives ∂ x i n applied at Similar results can be found for the x quantities substituting Equation (21) into Equation (12): Now we evaluate σ rel R * varying the x i values across the range reported in Equation (2) and by calculating E[n j ]. We see that σ rel R * is never higher than 10 −11 : In addition to the above uncertainties, the National Institute of Standards and Technology (NIST) estimates the error using the empirical Ciddor formula, that is unrelated to the uncertainties of the input atmospheric parameters. This takes the form of an adjustment to the formula accounting for various effects (please see https://emtoolbox.nist.gov/Wavelength/Documentation.asp#AppendixAV for detailed explanation and formula). By direct evaluation we estimate the error to be ∼10 −8 , while assuming higher values only at the very extreme values of the formula's validity.

General Case
In the general case, Equation (17) is not valid anymore and we need to consider the spatial distribution of n: R n(r)dr → Trajectory≡ L n(x(l))dl (28) and we can rewrite the equations above this way

Linear Variation of Atmospheric Parameters
In order to have a numeric estimation of the error σ rel R * in case of atmospheric parameters x that change values between the observer (the laboratory where the laser is positioned) and the target, we have chosen a specific value of the range R, a specific linear variation of x with light trajectory between the lab and the target, and performed numerically the integral in Equation (16). This provides an estimation of the real observed beat frequency f j .
Then we used Equation (12) to measure the range as per the multi-frequency method above exposed, and compared this prediction R * with the known value R. To this aim, we assigned to x (0) the x values as measured in the lab.
Similarly we used f j also to estimate the range in the single-frequency classical case R * * , using Equation (1), and compared it with the known R. For this case the refractive index used is calculated from the x in the lab conditions as well.
The results of the relative error are reported in Table 1. In the table we also report the example wavelengths of a Nd:YAG laser and its harmonics generated in nonlinear crystals. They produce beat frequencies of f beat = 13.5483, 13.5484, 13.5485, 13.5487, 13.5489 GHz for R = 10,101 m, with a chirp rate K = 201 MHz/µs, the same order of magnitude as in [27]. Similar relative errors are found for R = 100,301 m. Table 1 shows two sets of range measurements obtained from two linear variations of atmospheric parameters between lab and target. The two sets are chosen in order to have a large excursion in the lab-vs-target atmospheric value (case 1) and a more mild condition (case 2). The atmospheric parameters are shown in Table 2.
For λ = 1.0645 µm and more extreme variation (case 1), the refractive index between the lab and the target varies from 1.00020745 to 1.000342853. P s = 101,325 Pa is the standard atmospheric pressure. For completeness, we calculated the expected error using multiple wavelength generated by electro-optic modulator (see Appendix B). In this case the vicinity of the generated wavelengths limits the performance of the multi-wavelength approach. Furthermore also the error due to the numerical inversion of matrix C needs to be addressed. In fact matrix C is close to a singular matrix, due to the terms with the derivatives of n being very small compared with the terms with n. Hence the system is solved numerically, with a certain error on the solution. As described in Appendix A, this error turns out to be small compared with the errors above, and can be neglected. Table 1. Relative errors obtained from one ((R * * − R)/R) and multi-wavelength ((R * − R)/R) approach in two different atmospheric conditions (Case 1 and Case 2 as reported in Table 2). The first column of the table reports the number of wavelengths used by the new method, while the second column reports their values.
Multi-wavelength approach is widely used in interferometry and ranging fields with pulsed laser and recently has been applied to FMCW. In this context is fundamental a careful study of the refractive index variation due to atmospheric parameters (i.e., pressure, temperature, humidity and CO 2 concentration). In our approach we focused on the accuracy of length estimation and we performed a numerical simulation in order to compare with the single wavelength approach. The result is an improvement on the accuracy limited only by accuracy of Ciddor formula. In addition, a careful modeling of the atmospheric influence on a continuous wave laser beam (amplitude or frequency modulated or even not modulated) enables new applications in the communication, meteorological and metrological fields.

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

Appendix A. Errors Due to Numerical Evaluation
The system (11) is more precisely solved to find P by means of numerical methods than by calculating the inverse C −1 , due to the high degree of singularity of matrix C. There are different techniques to do so. The technique, among those we tried, that finds the solution P * of the system (11) producing the smallest error on the vector F is the LU decomposition method. If is the error on F produced by the decomposition we can write: To perform LU decomposition we used the numpy linalg.solve function that relies on the fortran LAPACK library-one of the standards for high precision and performance scientific computing [42,43], that returns the error ∼ 10 −12 m for every component. We can use to propagate the relative error R on the range R. In fact we can manipulate the first equation of system (11) as follows since R ∼ 10 4 m and (c 10 + ∑ l i=1 c 1i y i ) −1 ∼ 1. The method proposed in this paper has been implemented as a python script, with the usage of standard routines of the SciPy ecosystem [43]. No proprietary and specifically dedicated software was needed, neither any particular setting. The numerical evaluation of integral in Equation (16) has been interpolated over 1000 points implementing a classical trapezoidal algorithm. No significant change is observed for higher number of points. Scripts and materials can be made available upon request to the authors.

Appendix B. Multiple Frequencies Generation by Means of Electro-Optic Modulator
The use of an electro-optic modulator is an alternative approach with respect to a non linear crystal for the generation of multiple frequencies from one laser frequency. These frequencies are symmetrically spaced by a fixed gap around a central frequency. A gap value of 40 GHz is now available in commercial devices . This produces very close values of the wavelength and increases the relative error σ rel R * . In Figure A1 we report how the relative error decreases with the increase of the frequency gap for a central wavelength of 1.55 µm. We can see that only for high not experimentally feasible gaps the error is sufficiently reduced. A slightly smaller, but quite similar, error is seen using a wavelength of 1.0645 µm. Figure A1. σ rel R * relative error (absolute value) vs. frequency gap via electro-optic modulator (central wavelength 1.55 µm). Different colors correspond to different number of frequencies used by the proposed approach to range estimation. It can be seen that for 3 and 4 (red and green colors) frequencies the method reaches good levels of sensitivity only for unfeasible high gap values, while for 2 frequencies only (blue color) the relative error never significantly declines in the interval considered, looking like a constant value in comparison to the other cases.