1. INTRODUCTION

Lidar remote sensing could be an essential tool for future security and surveillance systems, especially for standoff detection and quantification of vapor phase chemical warfare agents (CWA) and toxic industrial chemicals (TICs) in the atmosphere. As most CWAs and TICs mainly absorb light in the extended long-wave infrared (LWIR) range between 6 and 14 µm [1], their long-range detection requires frequency-agile LWIR lidar systems. Such systems have already been reported in the past, generally relying on ${\rm CO}_2$ lasers [2–7]. However, the tunability range of ${\rm CO}_2$ lasers is restricted to discrete lines in the 9.2–9.8 and 10.1–10.8 µm bands, which does not allow addressing all relevant species for security applications. To overcome this limitation, the French Aerospace Lab (ONERA) has driven during these past years the development of a new lidar system, called mid-infrared chemical lidar (MICLID), widely tunable in the LWIR domain. The system is based on an all-solid-state pulsed laser source, featuring down-conversion parametric optics and amplifiers to reach the LWIR band. This paper presents the theoretical evaluation of the system’s performance for detection and quantification of several CWAs of interest and addresses several issues presented hereafter. The design and field-test performances of the built lidar system are the topic of the companion paper [8].

As a general tool for sensing a wide variety of chemical agents, MICLID has been designed to cover a large wavelength range between 7.5 and 10.7 µm. In a surveillance scenario, where the lidar has to evaluate potential threats from a given area, the lidar is required to cover the entire tunability range in a limited time, which implies using a limited set of probing wavelengths. In the case of MICLID, it was decided to build a 17-wavelengths set, distributed over the tunability range, to be covered within a 2 min alert time (${\sim}{7}\;{\rm s}$ per wavelength). This defines what is called hereafter the system’s “alert-mode” wavelength set. All of the wavelengths of the alert mode must be carefully chosen. They must coincide with absorbing lines of the chemicals under surveillance, but they also must avoid, or minimize, the interference level with natural atmospheric species. These issues are detailed in Section 2. In Sections 3 and 4, a performance model is built and used to evaluate the lidar precision for quantification of several species of interest. The analysis especially develops the case of an integrated-path differential absorption (IPDA) lidar mode, in which it is assumed that a topographic hard target backscatters the laser light at some distance behind the gas plume. Indeed, because CWA molecules are heavy, associated gas plumes propagate near the ground, such that any laser line of sight should be terminated by some hard target. In that case, only the path-integrated content of the chemical species is considered, which is enough to identify a threat. Cramer–Rao bounds (CRBs) for quantification of the species-integrated contents are used to evaluate the system’s performance. In Section 5, a numerical method is proposed in order to determine the optimal shortest set of wavelengths maximizing the lidar performance in a given time and for a given objective. The method is general and applicable to a variety of scenarios, and it is exemplified through a four-species optimization problem. Finally, in Section 6, the lidar receiver operating characteristics (ROC) curves are computed for the optimal set of wavelengths found in Section 5. These curves allow for evaluating the system’s detection performance, expressed in probability of detection (PD) as a function of the probability of false alarm (PFA) for a given detection scenario.

2. ATMOSPHERIC TRANSMISSION

As shown in the literature [1], a large number of CWA and TICs exhibit absorption features in the 6–14 µm band (so-called molecular fingerprint region). Below 7.5 µm and beyond 12 µm, it is known however that very strong absorption by water vapor generally precludes long-range detection. Therefore, MICLID was designed to cover the 7.5–10.7 µm range, which allows for addressing a large number of relevant gas species. In this section, we first illustrate how the atmospheric transmission puts stringent conditions when choosing wavelengths for a widely tunable surveillance lidar system in the LWIR range. We then present absorption transmission spectra of four specific CWA species that will be referred to in the following sections, illustrating how the alert-mode wavelength set coincides with those species absorption lines.

A. Transmission of Natural Atmospheric Species

As mentioned in the introduction, the lidar performance will be assessed assuming an IPDA mode, i.e., where a hard target is located at a distance $ Z_T $ from the lidar and backscatters a fraction of light toward the lidar receiver (Fig. 1).

 

Fig. 1. Measurement principle of the lidar operating in the IPDA mode. Two detection channels are used at each laser shot. The energy reference channel $ Q $-ref measures the shot-to-shot laser energy fluctuations (time-integrating the emitted pulse waveform), while the lidar signal channel $ Q $-lid measures the lidar signal (time-integrating the lidar echo waveform).

Download Full Size | PDF

 

Fig. 2. Two-way atmospheric transmission for a hard target at 500 m (${\rm H}_2{\rm O}$ 2%, ${{\rm CO} _2}$ 400 ppm, ${\rm CH}_4$ 1.8 ppm, ${\rm N}_2{\rm O}$ 0.32 ppm, ${\rm O}_3$ 0.027 ppm). Pink dots indicate the alert-mode wavelength set.

Download Full Size | PDF

 

Fig. 3. Two-way atmospheric transmission contributions for a hard target at 500 m (${\rm H}_2{\rm O}$ 2%, ${{\rm CO} _2}$ 400 ppm, ${\rm CH}_4$ 1.8 ppm, ${\rm N}_2{\rm O}$ 0.32 ppm, ${\rm O}_3$ 0.027 ppm).

Download Full Size | PDF

Figure 2 shows the total atmospheric transmission calculated for a hard target located at ${{ Z }_T} = 500\;{\rm m}$ (the transmission corresponds to a two-way path, totaling 1 km), while Fig. 3 shows the contribution (in transmission) of the main natural absorbing species in the 7.5–11 µm band. Those species, and their corresponding volume mixing ratios (VMR) chosen for the computation (corresponding to standard atmosphere composition), are ${\rm H}_2{\rm O}$ (${\rm VMR} = 2\%$ for a relative humidity ${\rm RH} = 71\%$ at temperature ${T} = 296\;{\rm K}$), carbon dioxide (${\rm CO}_2$, ${\rm VMR} = 400\;{\rm ppm}$), methane (${\rm CH}_4$, 1.8 ppm), nitrous oxide (${\rm N}_2{\rm O}$, 0.32 ppm), and ozone (${\rm O}_3$, 0.027 ppm). Small pink dots represent the alert-mode set of wavelengths. All of those plots have been computed with a line-by-line code relying on the high-resolution transmission (HITRAN) database [9]. For each absorption transition, the line-by-line calculation assumes a Voigt line shape truncated at ${25}\;{{\rm cm}^{- 1}}$ from the line center frequency. Consequently, the calculation accounts for cumulative wing effects only in the limit of the ${25}\;{{\rm cm}^{- 1}}$ truncation. In order to account for far-wing cumulative absorption effects beyond ${25}\;{{\rm cm}^{- 1}}$, specific models can be used, known as continuum absorption models. Among them, the Mlawer, Tobin, Clough, Kneizys, and Davis (MT_CKD) model [10,11] is an acknowledged reference. Light absorption by the MT_CKD continuum (version 2.5.2 released in 2011) is thus also shown in the ${\rm H}_2{\rm O}$ part of Fig. 3 and is commented on below. Light extinction due to Rayleigh scattering by molecules is neglected in this work, since the ${\lambda ^{- 4}}$ spectral dependency of Rayleigh scattering cross sections leads to very low values in the LWIR. Extinction due to Mie scattering aerosols is neglected for the computation with ${{ Z }_T} = 500\;{\rm m}$, but will be commented later when considering the case of a far-range lidar with ${{ Z }_T} = 5\;{\rm km}$ (see Fig. 4).

 

Fig. 4. Atmospheric transmission for a hard target at 5 km.

Download Full Size | PDF

Even though the LWIR domain is known to be a good atmospheric transparency window, it can be seen on Figs. 2 and 3 that even for a medium-range hard-target distance (500 m), strong ${\rm H}_2{\rm O}$ lines are ubiquitous. For a standoff lidar dedicated to CWA detection, it is thus mandatory to emit wavelengths that do not coincide with ${\rm H}_2{\rm O}$ lines. As shown on the figures, the 17 wavelengths chosen to build the “alert mode” comply with this requirement. They are also approximately equally distributed all over the 7.7–10.7 µm interval (with small irregularities to favor detection of specific CWA absorption features). The continuum absorption leads to non-negligible 5–10% absorption above 8 µm, with a rapid absorption increase below 8 µm. It is thus essential to take into account this continuum effect for a correct evaluation of the atmospheric transmission, hence, of the lidar performance. Figure 3 also shows that ${\rm CO}_2$, ${\rm CH}_4$, ${\rm N}_2{\rm O}$, and, to a lower extent, ${\rm O}_3$ can induce substantial absorption in some spectral bands. Attention must therefore be paid to these species when choosing the alert-mode wavelength set. ${\rm CO}_2$ absorption bands cover the intervals 9.2–9.8 µm and 10.1–10.8 µm, and, for a hard target at 500 m, peak at 15% around 9.3 µm. The alert-mode wavelengths are chosen to minimize ${\rm CO}_2$ absorption to 1–2% (residual absorption between the lines, due to cumulative wing effects). ${\rm CH}_4$ and ${\rm N}_2{\rm O}$ both exhibit very strong absorption features below 8.2 µm, and it is thus mandatory to choose wavelengths minimizing the interference level with these two species. For ${\rm N}_2{\rm O}$, it can be seen that the chosen wavelengths comply with this objective, though the residual absorption between ${\rm N}_2{\rm O}$ lines remains noticeable (7–8% here). The chosen wavelengths do not exactly minimize ${\rm CH}_4$ absorption, but a compromise must be made between ${\rm CH}_4$, ${\rm H}_2{\rm O}$, and ${\rm N}_2{\rm O}$ absorption lines, and Fig. 2 shows that the chosen wavelengths below 8 µm are optimal when considering the three species together. Finally, ${\rm O}_3$ can induce a small absorption in the 9.4–10 µm band (peaking at 4% at 9.5 µm). Here, again, the chosen wavelengths minimize the interference level to 1–2%.

Obviously, though the 17 wavelengths have been chosen so as to maximize atmospheric transmission, absorption by natural gas species will get worse if the lidar range is increased. Figure 4 shows the expected atmospheric transmission for a hard target at ${{ Z }_T} = 5\;{\rm km}$. It can be seen that the ${\rm H}_2{\rm O}$ continuum induces more than 40% absorption all over the LWIR band, while wavelengths under 8 µm are almost completely absorbed. Therefore, the water vapor absorption continuum is expected to have a strong impact on LWIR lidar performance and wavelength optimization strategies, especially at the far range and for gas species signing below 8 µm. For this far-range simulation, light extinction by aerosols scattering is also taken into account. Aerosols are assumed to be spherical and to belong to the “continental polluted” class defined in the online public access catalogue (OPAC) database [12], with a mass concentration of ${50}\;\unicode {x00B5}{\rm g}/{\rm m}^3$. It can be seen that the aerosol extinction remains limited (10% extinction at peak) compared to the water vapor absorption continuum for instance.

In conclusion, choosing a good wavelength set for standoff detection of CWA on a wide spectral range in the LWIR domain requires careful attention to the atmospheric transparency properties, as many wavelengths could be strongly absorbed by natural species. In principle, absorption by those molecules could be anticipated and pre-corrected through ancillary in-situ measurements or by regularly calibrating the lidar response in the absence of target gas. In such conditions, absorption by the natural atmosphere would not bring any bias on target gas measurements, yet they will affect the measurement precision by impacting the signal-to-noise ratio (SNR) recorded at each wavelength.

B. Transmission of Chemical Agents

Though the MICLID lidar has been designed as a general tool for CWAs and TICs detection in the LWIR band, we will focus, from here and for the following sections, on a measurement scenario involving four specific target species. Three of them are important CWAs: sulphur mustard (NATO-designated by letters HD), Sarin (GB), and V-base nerve agent (VX). The fourth species is a simulant often used in infrared spectroscopy, called triethyl phosphate (TEP). Absorption cross-section spectra of GB, HD, VX, and TEP can be found in the literature [1,13]. From these data, it is possible to compute the two-way transmission of a gas mixture containing $ S $ different species using the Beer–Lambert law:

 

Fig. 5. Transmission spectra of GB, HD, VX, and TEP for their respective TSLs. Green circles represent the optimized set of five wavelengths selected by the stepwise wavelengths reduction algorithm, and blue numbers in the bottom figure indicate the iteration loop number at which the other wavelengths are eliminated (see Section 5).

Download Full Size | PDF

3. PERFORMANCE MODEL

The IPDA technique consists of emitting a series of laser pulses for each wavelength of the alert-mode set and in recording and time-integrating the signals backscattered by a hard target at distance $ Z_T $, as shown in Fig. 1. Computing, for each pulse, the ratio between the integrated signal and the emitted pulse energy, gives access to a signal proportional to the atmospheric two-way transmission between the lidar and the hard target. Pulse averaging is generally used to improve the SNR before stepping to the next wavelength. At the end of the allocated measurement time, when all of the desired wavelengths have been transmitted and the lidar echoes collected, a multi-wavelength estimation technique can finally be employed to quantify each compound under surveillance. This section first establishes the expression of the measurement variance of the recorded signal for a single wavelength. Then, it presents the multi-wavelength equations that can be used to calculate the expectable variance on quantified CL products for the species under surveillance. The treatment borrows from previous works in the field [14–17] and is summarized here for self-consistency of this work.

A. Signal Variance Model for a Single Wavelength

As illustrated in Fig. 1, for each emitted pulse, the reference and lidar receivers deliver photocurrent waveforms ${i_{{\rm ref}}}(t)$ and ${i_{{\rm lid}}}(t)$, centered at time ${t_{{\rm ref}}} = 0$ and ${t_{{\rm lid}}} = {{2Z}_T}/c$, respectively, where $ c $ is the speed of light and $ Z_T $ the hard-target distance. Time-integration of these signals over durations $ T_{\rm ref} $ and $ T_{\rm lid} $ comprising whole waveforms yields the following signals (in number of electrons):

The random variable $ N $ accounts for shot-to-shot laser energy variations. Fluctuations greater than 10% or 20% (sometimes much more) are common in high-energy pulsed laser sources. For this reason, the emitted shot-to-shot energy is measured using the reference channel. The speckle modulation index ${\nu _s}$ [14,15] describes the energy fluctuations of the light collected on the telescope aperture, as a result of the varying interference pattern generated by laser scattering on the hard target. This pattern is sensitive to the hard-target roughness at the wavelength scale. It can therefore vary very quickly, for instance, if the hard target is moving, shaking, or in case of beam pointing jitter. Speckle noise is not accounted for in the reference channel, because in this channel the laser beam spatial coherence can in principle be conserved (no interference pattern on the reference detector). The baseline modulation index ${\nu _B}$ reflects possible time fluctuations of the $ B_0 $ value during the system’s measurement time (a total of 2 min for 17 wavelengths). For instance, such fluctuations could come from hard-target reflectivity variations ($ \kappa $ coefficient), in case of hard-target motion, laser-induced or turbulence-induced beam pointing fluctuations, or from atmospheric transmission variations (${T_{{\rm ATM}}}$ coefficient) like time-varying contents of water vapor or aerosols. Finally, $ n $ and $ m $ represent the lidar channel and reference channel detector noises.

For the following, let us note $\langle X \rangle$, $ \sigma_X $, and ${{\rm SNR}_X} = \langle X \rangle /{\sigma _X}$ as the statistical mean, standard deviation, and SNR of a random variable $ X $. The number of emitted photons per pulse $ N $ has a mean $\langle N \rangle$ and variance $\sigma _N^2 = \langle {{N^2}} \rangle - {\langle N \rangle ^2}$. Speckle and baseline modulation indices are unitary random variables with mean $\langle \nu \rangle = 1$ and variances $\sigma _\nu ^2 = \langle {{\nu ^2}} \rangle - 1$. Detection noises $ n $ and $ m $ are zero-mean random variables with variances $\sigma _n^2$ and $\sigma _m^2$. In the ${Q_{{\rm lid}}}$ channel, $ N $, ${\nu _S}$, ${\nu _B}$, and $ n $ are statistically independent one from each other since there are no physical correlations between these fluctuations sources. For the same reason, $ N $ and $ m $ are also independent in the $ Q_{\rm ref} $ channel. Computing the signal ${R} = {\rm ln}({{Q}_{{\rm lid}}}/\!{{Q}_{{\rm ref}}})$ allows both normalizing the laser shot-to-shot fluctuations and linearizing the Beer–Lambert equation Eq. (1). The mean and variance of $ R $ can be calculated using a Taylor expansion of the logarithm of a ratio of random variables. We obtain [18]:

At the numerator, we have $\langle {{Q_{{\rm lid}}}} \rangle = N \cdot {B_0} \cdot {T_{{\rm CWA}}}$ in the lidar channel and $\langle {{Q_{{\rm ref}}}} \rangle = N\alpha$ in the reference channel. At the denominator, the first term corresponds to the signal shot noise, with $ F $ the excess noise factor possibly induced by the detection chain (for instance, the Fano factor of an avalanche photodiode and/or other amplifiers), the second term corresponds to the background signal shot noise, and the third term corresponds to the electronic circuit noise. In each channel, the background signal integrated during time $ T $ can be modeled by

The speckle term ${\rm SNR}_{{\nu _S}}^2$ is equal to the number of spatio-temporal speckle cells averaged on the detector surface (or, equivalently, on the telescope aperture) during the waveform integration time. This number can be approached [15] by

The temporal speckle averaging number $ m_T $ is the number of independent speckle patterns created by the laser or lidar beam on the detector during the integration time $ T $. It is the ratio between the incident waveform duration and its coherence time ${m_T}\approx {t_{{\rm wav}}}/{t_{{\rm coh}}}$. Emitted pulses exhibiting Fourier-transform-limited spectral contents typically have coherence times equal to the waveform duration, leading to ${m_T}\approx 1$. On the contrary, chirped pulses have lower coherence times, leading to higher $ m_T $ values.

The baseline fluctuations level given by ${\rm SNR}{\nu _B}$ is not easy to anticipate or to model analytically. Indeed, it is expected to be very dependent of experimental conditions such as the total measurement duration, the hard-target motion, atmospheric dynamics (humidity, aerosols…), beam pointing fluctuations, etc. A sensitivity study is conducted in Section 4 to evaluate its potential impact.

Finally, pulse averaging can be performed to reduce the $ R $-signal variance. In principle, the pulse averaging number $ N_{\rm avg} $ could be different for each emitted wavelength, but, for simplicity, it will be considered constant for all wavelengths throughout this paper. The averaging effect can be simply modeled by multiplying the values of ${{\rm CNR}_{{\rm Qlid}}}$, ${{\rm CNR}_{{\rm Qref}}}$, ${\rm SNR}{\nu _S}$, and ${\rm SNR}{\nu _B}$ by ${({N_{{\rm avg}}})^{1/2}}$. Doing so is physically justified for the CNR terms, because they represent detection noises that are statistically independent from pulse to pulse. On the contrary, pulse-to-pulse independence of the speckle and baseline noises may not be strictly verified. For instance, if the hard target and the incident laser beam footprint were perfectly motionless, the speckle pattern should remain identical from pulse to pulse (acting as a bias rather than a noise in that case). But in practice, target micro-motions at the wavelength scale and beam pointing instabilities (laser pointing jitter, turbulence-induced beam wandering…) should often combine to make the speckle pattern a random noise with short time correlation, especially for far-range lidar systems. For simplicity, though partial speckle time correlation may subsist from pulse-to-pulse in some real situations, complete decorrelation will be considered in this paper. The same approximation will be made for the baseline noise, since it looks reasonable to consider again that pulse-to-pulse beam pointing fluctuations could induce complete time decorrelation of the baseline noise induced by hard-target reflectivity fluctuations. However, a baseline noise component caused by atmospheric dynamic fluctuations (water vapor or aerosols content fluctuations) would certainly lead to partial time correlation from pulse to pulse. This effect could be assessed in future dedicated studies by introducing coherence time parameters.

B. Estimate Variance for Gas Species CL Products

The lidar system is designed to emit sequences of ${M} = 17$ successive wavelengths within 2 min in the alert mode, and, for each wavelength, the analytic model developed in the preceding section applies. In matrix notation, Eq. (4) can be written as

Quite often, the baseline is a smooth function of wavelength that can be modeled by a weighted linear function. For instance, a $ P $-order polynomial can be assumed, such that ${K_i} = \sum_{k = 0}^P {{a_k}\lambda _i^k}$, where all of the $ a_k $ coefficients are supplementary unknown parameters. In the following, we note $\boldsymbol{\theta} = (C{L_1}, \ldots C{L_S},{a_0} \ldots {a_P})$ as the vector containing all of the unknown parameters. In best cases, the mean $ \boldsymbol{K} $ baseline can be measured through a calibration measurement (sometimes up to a scaling constant such that only $ a_{0} $ has to be considered as a supplementary unknown parameter).

In Eq. (12), $ \boldsymbol{R }$ is a multivariate random variable, whose covariance matrix is noted as ${\boldsymbol{\Gamma} _{\boldsymbol{R}}}$. The diagonal elements of ${\boldsymbol{\Gamma} _{\boldsymbol{R}}}$ are given by Eq. (5) multiplied by $N_{{\rm avg}}^{- 1}$. Non-diagonal elements could be different from zero if speckle and baseline fluctuations were correlated between the wavelengths successively emitted by the lidar, but as explained in the preceding section, only the complete time-decorrelation case is treated in this paper, and consequently ${\boldsymbol{\Gamma} _{\boldsymbol{R}}}$ is assumed to be diagonal. Assuming also that each random variable in $ \boldsymbol{R} $ follows a normal distribution (see the discussion in Section 6 on that point), then the maximal likelihood (ML) estimator of $ \boldsymbol{\theta }$ and its associated covariance matrix ${\boldsymbol{\Lambda} _{\boldsymbol{\theta}}}$, are given by [14]

The estimate variance for the parameter ${\theta _j}$ can then be computed from the $ j $th diagonal element of ${\boldsymbol{\Lambda} _{\boldsymbol{\theta}}}$:

It can be shown that ${\boldsymbol{\Lambda} _{\boldsymbol{\theta}}}$ has the same expression as the inverse of the system’s Fisher information matrix (FIM) [20]. Thus the diagonal of ${\boldsymbol{\Lambda} _{\boldsymbol{\theta}}}$ identifies with the system’s CRB, representing minimal variances reachable by an unbiased estimator. In the following, we shall only focus on the CRBs applying to the CL parameters (${\theta _1}$ to ${\theta _4}$). In order to assess the system’s performance compared to observational objectives, we also introduce the neasurement quality factor (MQF) for each species by

Table 1. Laser and Lidar Parameters

View Table

4. SIMULATION OF PERFORMANCES

A. Lidar Parameters

The lidar parameters for performance simulations are given on Table 1 and are based on realistic design considerations. Most of them are also supported by laboratory measurements, as discussed in the companion paper [8]. The laser is assumed to be a tunable single-mode optical parametric oscillator (OPO)-based pulsed transmitter. Using a nested-cavity OPO scheme, it is possible to obtain single longitudinal mode (SLM) emission without seeding the laser cavity. The laser linewidth is considered Fourier transform limited (${\lt}{100}\;{\rm MHz}$) by the pulse duration (10 ns). It can thus be neglected compared to the broad cross-section lines of the target gas compounds [several tens of gigahertz (GHz)]. This assumption also implies a unit temporal speckle number (no frequency chirp, hence, ${m_T} = 1$ in Eq. (10)). The laser half-divergence angle and the beam waist radius on the emission optics are set to 1 mrad and 6.2 mm at 10 µm. This corresponds to a beam quality factor of ${{ M }^2} = 2$. Beam waist and ${M^2}$ values are assumed to be independent from the wavelength, as discussed in Eq. (11). With such parameters, the average number of spatial speckle cells is about ${\rm Ms}\approx {260}$; hence, the speckle ${{\rm SNR}_{\nu S}}$ is around 16. A telescope diameter of 10 cm, featuring a central obstruction mirror of 2 cm in diameter (used for on-axis beam emission), is also assumed, and the two-way optical losses are supposed to reach 50% due to imperfect coatings of emission and reception optics or other imperfections in light collection efficiency. The hard-target albedo is taken to 5% (typically corresponding to a black piece of fabric). The overlap factor is 100% at the target distance (far-field focused lidar), and the target depth profile is assumed to reach 5 m here (typically a multi-plane target like trees). The detector is assumed to be a 1 mm active diameter HgCdTe photodiode, integrating the incoming light over a 1–11.5 µm spectral broadband, with 50% quantum efficiency (sensitivity of 4 A/W at 10 µm), 60° field of view, and a NEP of ${1.1}\;{{{\rm pW}/{\rm Hz}}^{1/2}}$ (or ${D}^* = 8.3^{*}{{10}^{10}}\;{\rm cm} \cdot {{\rm Hz}^{1/2}}/{\rm W}$). An excess noise factor of ${F} = 1.5$ and a bandwidth of 5 MHz are assumed for the detection and amplifier chain. With those parameters, the electrical waveform delivered by the ${Q_{{\rm lid}}}$ channel for each emitted laser pulse has a temporal FWHM of 75 ns (mainly limited by the detector bandwidth). The waveform integration time $ T $ is thus chosen to be 190 ns; hence, the integration bandwidth $ B_T $ is 2.7 MHz. It is assumed first that the full 17-wavelength alert-mode sequence is used. With the alert time being set to 2 min, the accumulation time per wavelength is 7 s. The baseline noise is first neglected (infinite ${{\rm SNR}_{\nu B}}$). The reference channel CNR is assumed to be very good (${{\rm CNR}_{{\rm ref}}} = 100$). Finally, as for the laser pulse energy and repetition rate, we evaluate two different power distribution strategies. Assuming an average laser power varying in the range of 1–100 mW, we refer to it as a “type-1 laser,” a laser with a fixed repetition rate (50 Hz) and growing pulse energy (from 20 µJ to 2 mJ), and as a “type-2 laser,” a laser with a fixed pulse energy (20 µJ) and growing pulse repetition frequency (from 50 to 5000 Hz).

 

Fig. 6. MQFs for GB, HD, VX, and TEP, as a function of the laser average power and laser type (type 1/type 2) for a hard-target distance of 800 m (left) or 1600 m (right).

Download Full Size | PDF

B. Lidar Performance

Based on the lidar parameters listed in Table 1 and model equations presented in Section 3, Fig. 6 shows the computed MQFs for GB, HD, VX, and TEP as a function of the laser average power and the laser type (type 1/type 2), for a hard-target distance of 800 m (left) and 1600 m (right). The situation at 800 m illustrates a classical tradeoff in lidar design. For the type-1 laser (50 Hz and growing pulse energy), the MQFs increase quickly with growing pulse energy until some saturation occurs. The saturation is due to the speckle noise. Indeed, at some point, the pulse energy is high enough for the detection noises (shot noise and detector noises) to be negligible compared to the speckle noise. Since the latter is of a multiplicative nature, it is always proportional to the pulse energy. Therefore, once the speckle ${{\rm SNR}_{\nu S}}$ becomes the limiting value in Eq. (5), the performance cannot be improved further by increasing the pulse energy. This is the reason why the type-2 laser performs better at higher laser power for a fixed distance. Indeed, the type-2 laser has a fixed pulsed energy, and a growing repetition rate. In that case, the detection performance improves by a speckle averaging process. The gain in SNR is proportional to the square root of the number of pulses as long as the assumption (discussed in Section 3) of uncorrelated speckle patterns for successive pulses can be made (otherwise the gain may be slower). In general, and especially for short-range IPDA applications, attention must therefore be paid to the expected speckle level in order to decide which of the type-1 or type-2 laser design strategies is better suited to the measurement objectives. The left plot in Fig. 6 shows, for instance, that for GB measurement at 800 m, 50 mW lasers with 20µJ/2500 Hz or 1 mJ/50 Hz would have similar performances. Beyond that point, the performance only improves with power with a type-2 laser. However, such tradeoff is highly dependent from the target distance. The right plot shows that for longer-range systems, type-1 lasers perform better than type 2. As long as the shot-noise and detector noise limit the performance, it is worthwhile to increase the laser pulse energy instead of the repetition rate.

In the following, we will consider only type-1 lasers in order to evaluate performances for long-range scenarios, since most of the previously reported CWA lidar systems use lasers of this kind (for instance 1 mJ at 20 Hz in [21]). The MQFs obtained show that at 1600 m, the target species should be measured with uncertainties at least 10 times below their TSL for a laser power above 30 mW for a 50 Hz, 0.6 mJ laser system. The easiest-to-detect target gas is GB (MQF of 30 for 30 mW on the right plot), while the most difficult is TEP (MQF of 10 for 30 mW).

 

Fig. 7. MQFs for GB, HD, VX, and TEP, as a function of the baseline SNR, for a hard-target distance of 1600 m and a 50 mW type-1 laser (1 mJ/50 Hz).

Download Full Size | PDF

Figure 7 shows the impact of possible baseline fluctuations for a hard target at 1600 m and a 1 mJ/50 Hz laser. As discussed in Section 3, the baseline might fluctuate in practice during the wavelength scan performed by the lidar because of variations of the water vapor content in the line-of-sight, the hard-target albedo, or the field overlapping function. In this example, a baseline ${{\rm SNR}_{\nu B}}$ limited to 10 makes the MQF drop by approximately 40% compared to an ideal non-fluctuating baseline. For precise performance assessment of the instrument, it will thus be important in future work to evaluate experimentally the baseline ${{\rm SNR}_{\nu B}}$. A way to avoid baseline fluctuations would be to use series of wavelength pairs, each pair being composed of two close wavelengths serving as local ON/OFF differential absorption lidar (DIAL) pairs. This approach (suggested also in [14]) would allow correlating the baseline noise in each wavelength pair and suppress it by computing their signal ratio. However, the spectral distance within the pair should be sufficient to distinguish some differential absorption feature of the target gas as well, and, since CWAs exhibit very broad absorption lines, this could require a non-negligible spectral distance that the laser would have to cover quickly. Note also that the baseline issue is expected to be more critical when considering an IPDA measurement mode compared to a range-resolved DIAL mode. Indeed, in the former case, the complete atmospheric column cumulates in the signal, and consequently, time variations of the cumulated water vapor content are more likely to occur and to fluctuate than in a limited-range gate. This issue could, therefore, be alleviated in future range-resolved systems.

 

Fig. 8. Evolution of the MQF of each CWA as the wavelength sequence is reduced. In the left plot (objective A), the total measurement time reduces in accordance with the number of wavelengths. In the right plot (right), the total measurement time is kept unchanged (2 min), hence the measurement time per wavelength increases accordingly with the number of suppressed wavelengths.

Download Full Size | PDF

5. WAVELENGTH SEQUENCE REDUCTION

In the above sections, the alert-mode scenario defines a set of 17 wavelengths in order to make a general cover of the LWIR domain. However, in many scenarios, only a limited number of CWAs would be expected. In such a case, some wavelengths of the alert-mode set might be useless, and the lidar could be optimized. In this section, we consider a fictional situation where only the four above-mentioned CWAs (GB, HD, VX, and TEP) have to be detected. Theoretically, if the background spectrum can be considered as known up to a single constant, the wavelength sequence could be reduced from 17 to 5. The five estimated parameters would be the CL content of the four potential agents plus the baseline level parameter $ a_{0} $. Reducing the wavelength sequence from 17 to 5 would potentially reduce the total measurement time from 2 min to 35 s, assuming equally distributed measurement time per wavelength. Another possible strategy would consist of increasing the measurement time per wavelength, while keeping the total measurement time unchanged in order to improve the [CL content measurement accuracy. Such optimization issues have been partly treated before by our team in a simpler case [22] (two-compound mixture and constrained wavelength choice). Here, we expand the discussion by considering a four-compound optimization problem and a lidar system able to emit any subset of wavelengths among those that build the general “alert-mode” set.

In order to determine which sequence of five wavelengths, within the alert-mode set, would be the best one, a specific stepwise algorithm has been developed. This algorithm aims at suppressing, at each step, the wavelength that is degraded the least and the smallest MQF factor among the four CWAs. Let us call objective “A” the objective of time reduction. The measurement time for each wavelength is $\Delta t = 7\;{\rm s}$ (350 pulses), and the total measurement time is ${T} = {N} \cdot \Delta t$ for an $ N $-wavelength sequence. Let us call objective “B” the objective of accuracy improvement at constant measurement time. In that case, the measurement time for each wavelength is $\Delta t = {120/\!N}\;{\rm s}$, and the total measurement time is constant at ${T} = 120\;{\rm s}$ (2 min) for an $ N $-wavelength sequence.

We ran the algorithm for a far-range case (hard target at 1600 m) so that wavelengths below 8 µm absorbed by the water vapor absorption continuum might not be retained as “good” wavelengths. A type-1 laser of 1 mJ/50 Hz and no baseline fluctuations are assumed. MQF evolutions with the wavelength sequence length are given on Fig. 8 for objective A (left plot) and objective B (right plot). All MQFs for the full 17-wavelength sequence are consistent with the values given in Fig. 6 (right) for a 50 mW type-1 laser. The smallest MQF for a full sequence is obtained for TEP, which implies that the wavelength suppression choice is governed by its effect on the MQF of TEP during the first iteration loops. After reducing from 17 to 7 wavelengths, following the suppression sequence labeled in Fig. 5 (bottom), VX becomes the smallest MQF among the four agents, therefore it drives the algorithm decision making about further wavelength suppression.

For objective A (time reduction), Fig. 8 (left) shows that reducing the sequence from 17 wavelengths to five wavelengths reduces the MQF by no more than 30% in the worst case (GB), while being faster by a factor 3.5. The measurement performances for VX and TEP are marginally affected by the wavelength suppression process, which indicates clearly that the remaining wavelengths contain all of the relevant information for those two species. Note that the two wavelengths under 8 µm are suppressed by the algorithm, since they contain only small absorption features for VX and TEP, and also because the water vapor continuum strongly attenuates them.

For objective B (accuracy improvement), Fig. 8 (right) shows that reducing the wavelength sequence while keeping the total measurement time unchanged is globally beneficial for the MQFs of all four CWAs (${+}{58}\%$ for TEP, ${+}{36}\%$ for HD, ${+}{27}\%$ for GB, ${+}{18}\%$ for VX). It also shows that for VX, the optimal wavelength sequence is obtained for six wavelengths. Further reducing the sequence to five wavelengths reduces the MQF for VX, though the suppressed wavelength is chosen so that the reduction is minimized. In the same way, for GB and HD, seven-wavelength sequences are optimal.

The five-wavelength sequence selected as the best one is shown on Fig. 5 (green circles). It is the same for objectives A and B. Indeed, only the number of pulses per wavelength distinguishes objectives A and B, which acts as a simple scalar factor of the covariance matrix ${\boldsymbol{\Gamma} _{\boldsymbol{R}}}$. Therefore, even though MQFs differ in absolute values during A and B iteration loops, their relative weights remain identical, leading to the same suppression choices.

6. RECEIVER OPERATING CHARACTERISTICS

In the preceding sections, we have assessed the system’s performance in terms of measurement precision of CWA-integrated contents. However, security lidar systems also have to operate in binary detection mode, i.e., provide a “yes or no” answer to the question “is this CWA present along the line of sight ?” Performance of this kind of detection systems is usually assessed by computing its ROC curves. Initially used for radar systems [23], a ROC curve typically plots the PD of an instrument as a function of its PFA. In this section, only the simplest case problem is detailed, where the system must decide whether a single compound is present in the line of sight or not. Multi-compound decision problems, involving multi-dimensional approaches such as those discussed in [14], are left out of the scope of this paper.

Let us suppose for instance that the system has to answer the question “is GB present in the line of sight or not?,” and let us call $ H_1 $ the “yes” hypothesis and $ H_0 $ the “no” hypothesis. To make a yes/no decision, we further assume that a classical likelihood ratio test (LRT) is applied to the measured $ \boldsymbol{R} $ vector Eq. (12). The LRT is defined by

It is the logarithm of the ratio of conditional probabilities to observe the actual signal $ \boldsymbol{R}$ under hypotheses $ H_1 $ and $ H_0 $. As a function of signal-dependent probabilities, the test value $ T $ is itself a random variable. If the test value is above (below) a certain threshold, we may want to declare that $ H_1 $ ($ H_0 $) is true, with more or less risk for a bad decision depending on the chosen threshold value. In order to evaluate the risk of providing false alarms (declare $ H_1 $ while $ H_0 $ is true) and the ability of the system to perform good detections (declare $ H_1 $ while $ H_1 $ is true), the conditional statistic of the LRT must be derived. Here, they are derived under simplifying assumptions that (1) all observable variables $ R_i $ are Gaussian random variables, (2) the standard deviation of a variable $ R_i $ is the same under $ H_1 $ and $ H_0 $, and (3) the baseline value and measurement standard deviations are both precisely known.

Assumption 1 should hold “near true” in our case. Noise affecting $ R_i $ values takes into account speckle noise, background shot noise, signal shot noise, detector noise, and “normalization noise” (noise added when dividing the received signal with ${Q_{{\rm ref}}}$). Detector noise is generally Gaussian, and speckle noise can be considered as Gaussian in direct detection because the number of speckle cells averaged over the telescope surface is very high (${\rm Ms} = 260$ in Section 4). Signal shot noise and background shot noise can also be considered Gaussian (instead of Poisson distributed) because the received average photon number is high enough (${\gt}\!{10}$) during the integration time. Dividing the signal by ${Q_{{\rm ref}}}$ will not affect the probability function much, as long as ${Q_{{\rm ref}}}$ is measured with a high SNR (an SNR of 100 is assumed here). Finally, computing the logarithm of a Gaussian variable should also change the probability function, but, in our case, because the signal SNR is quite high, this will remain a second-order effect. The resulting statistic of each $ R_i $ can therefore be considered as Gaussian in the first approximation. Assumption 2 is not true strictly speaking. Indeed, standard deviations for speckle noise and signal shot noise are both signal-dependent and should therefore be different under $ H_1 $ and $ H_0 $. However, this assumption can still be made because we consider low optical absorption, (${\lt}{15}\%$). Whether a CWA is present or not, it will not change the received signal very much (which is by the way the reason why a decision test is needed). Assumption 3 is also reasonable, since a gas surveillance lidar would be expected to measure the baseline alone most of the time (no gas threat) by performing repetitive wavelength scans. The baseline mean and standard deviation should thus be accurately known in many practical cases.

Given these assumptions, the LRT conditional probability density functions (pdf) under $ H_1 $ and $ H_0 $ can be obtained. If the goal of the decision test is just to decide whether the signal consists of the baseline alone ($ H_0 $) or the baseline plus a known CL value of the target gas (during a controlled detection test for instance), then the LRT conditional pdfs are given by simple Gaussian functions [20]. We obtain

 

Fig. 9. Conditional probabilities of the LRT tests under H0 and H1 for the detection of GB concentrated at TSL/10. Left: controlled test conditions (the target gas concentration is known); right: blind test conditions (the target gas concentration is estimated by an ML estimator).

Download Full Size | PDF

In Fig. 9 (left), we show the histograms of the conditional probabilities ${P}({{T}^\prime}|{H}_0)$ and ${P}({{T}^\prime}|{H}_1)$ that can be obtained in the case of a controlled GB detection test, at concentration TSL/10, for a hard target at 1600 m distance, with the optimal five-wavelength set identified in Section 5, and a total measurement time of 35 s (strategy A in Fig. 8). The distance $ d $ calculated from Eq. (21) in this case is $d = 4.45$. To build each histogram, 10,000 $ R $-vector random samples have been generated that verify the given covariance matrix ${\Gamma _R}$. Then, the $ T $ and ${{ T }^\prime}$ test values have been computed for each vector sample, assuming that the baseline and the gas quantity under test are known. The expected ${{ T }^\prime}$ test statistics [${N}({0},{1})$ and ${N}({d},{1})$] are correctly found.

However, in general, the CL value of the target gas to detect is not known from the system. In such blind-test operations, the goal of the decision test is to decide whether the signal consists of the baseline alone ($ H_0 $) or the baseline plus an unknown CL value that has to be estimated from the data at hand. Assuming that a ML estimator is used, such as given in Eq. (13), the LRT conditional pdf are given in [24]:

In this case, the conditional pdf under $ H_0 $ is a central chi-square distribution with one degree of freedom, while the conditional pdf under $ H_1 $ is a non-central chi-square distribution with one degree of freedom and with the decenter parameter $ d $ given again by Eq. (21). In Fig. 9 (right), we show the histograms that are obtained for a blind GB detection test at concentration TSL/10 (and the same other conditions as formerly stated). Again, the histogram shapes agree very well with the calculated pdfs Eqs. (22) and (23). We have also checked that the histogram for the ML estimates of the CL values has a regular Gaussian shape with a correct mean at TSL/10 and a standard deviation matching with the CRB calculated from Eqs. (13)–(16).

From the conditional pdf of the LRT under $ H_0 $ and $ H_1 $, it is possible to calculate the PFA and PD. For the controlled test case, where the pdf are Gaussian, these probabilities can be written explicitly as

Because the erfc* function is monotonic, it is equivalent to formulating a test threshold value $ \varepsilon $ or a PFA value. The PD can thus be expressed directly as a function of the PFA. Using the Matlab error function erf (slightly different from the erfc* definition given in Eq. (26)) and reverse function erfinv, the ROC can be calculated using the following equation:

Similarly, the PFA and PD values can be computed for the blind-test case by numerically integrating the statistics given by Eqs. (22) and (23) above a certain threshold.

The calculated ROC curves for the controlled detection test and the blind detection test are shown in Fig. 10. As can be expected, the ROC curve in blind test conditions is lower than for the controlled test, since not knowing the target gas CL value brings additional uncertainty. The curves are close to each other however, which can be explained by the fact that a high MQF is reached for GB (${\rm MQF} = 25$) in the assumed measurement conditions. Therefore, the cost for estimating the CL value is not very high. For a ${\rm PFA} = {{10}^{- 4}}$, we obtain, for instance, the PD at a 70–75% concentration level of TSL/10. Note that a PFA level of ${{10}^{- 4}}$ corresponds to an acceptable rate of about 0.1 false alarm per day, assuming that the lidar updates test decisions every 105 s (time needed to test GB, HD, and VX with 35 s per test, as assumed in this simulation). Such an update time is realistic and would enable a human being to take appropriate protection measures against a toxic cloud detected at 1 km and pushed by a 10 m/s wind.

 

Fig. 10. ROC curves for GB detection concentrated at TSL/10 in controlled test conditions and in blind test conditions.

Download Full Size | PDF

Similar ROC curves are obtained for HD at a concentration level of about TSL/6.5 and for VX at TSL/4.5. These results suggest that in the assumed scenario, good detection rates should be reached by the lidar system even for concentrations lower than the TSL. Ideally, such ROC curves would have to be checked experimentally in future work, for instance, by applying an LRT test to a large number of independent experimental data sets.

7. CONCLUSION

The expected performance of a multi-wavelength LWIR lidar has been assessed for IPDA measurement of CWAs. For long-range applications, it has been illustrated that high pulse energy and low repetition frequency are the best way to distribute the laser power. A detailed performance model for gas quantification has been developed, indicating that measurements with good sensitivity can be obtained at a long distance with reasonable laser power (below 100 mW) within 2 min. A wavelength-reduction algorithm has also been developed in order to determine the best subset of wavelengths (within a predefined larger set) for a given measurement objective. The presented method is general and can be applied to a variety of scenarios. Very significant performance improvements (in time duration or in measurement accuracy) can be obtained by optimizing the wavelength subset. ROC curve calculations have also been detailed for controlled (blind) detection conditions, where the gas concentration of the test is known (unknown). In each case, the anticipated statistics of LRT conditional probabilities agree very well with the histogram shapes that can be obtained from simulated signals. These ROC curves are useful for evaluating the capacity of a lidar instrument to be used as a future alert system.

In perspective, this theoretical work can be improved and completed in many ways. For instance, noise properties in the IPDA mode could be refined and confronted with real measurement noises. Time correlation parameters of speckle and baseline noises could be taken into account to improve the performance model accuracy when averaging pulses. ROC curve calculations may also be enriched by treating more complicated cases like multiple gas hypothesis tests or by considering as unknown the baseline properties (mean and standard deviation, see assumptions in Section 6). Finally, the proposed model could be straightforwardly extended to address the case of a range-resolved lidar, relying on aerosol backscattering instead of a hard target. In such a mode, the detection precision would be degraded, but the lidar could better localize a threat, and cumulative effects of the fluctuating atmosphere (baseline noise) could be alleviated.