Extending water vapor measurement capability of photon-limited differential absorption lidars through simultaneous denoising and inversion

The micropulse differential absorption lidar (MPD) was developed at Montana State University (MSU) and the National Center for Atmospheric Research (NCAR) to perform range-resolved water vapor (WV) measurements using low-power lasers and photon-counting detectors. The MPD has proven to produce accurate WV measurements up to 6 km altitude. However, the MPD's ability to produce accurate higher-altitude WV measurements is impeded by the current standard differential absorption lidar (DIAL) retrieval methods. These methods are built upon a fundamental methodology that algebraically solves for the WV using the MPD forward models and noisy observations, which exacerbates any random noise in the lidar observations.

The work in this paper introduces the adapted Poisson total variation (PTV) specifically for the MPD instrument. PTV was originally developed for a ground-based high spectral resolution lidar, and this paper reports on the adaptations that were required in order to apply PTV on MPD WV observations. The adapted PTV method, coined PTV-MPD, extends the maximum altitude of the MPD from 6 to 8 km and substantially increases the accuracy of the WV retrievals starting above 2 km. PTV-MPD achieves the improvement by simultaneously denoising the MPD noisy observations and inferring the WV by separating the random noise from the non-random WV.

An analysis with 130 radiosonde (RS) comparisons shows that the relative root-mean-square difference (RRMSE) of WV measurements between RS and PTV-MPD exceeds 100 % between 6 and 8 km, whereas the RRMSE between RS and the standard method exceeds 100 % near 3 km. In addition, we show that by employing PTV-MPD, the MPD is able to extend its useful range of WV estimates beyond that of the ARM Southern Great Plains Raman lidar (RRMSE exceeding 100 % between 3 and 4 km); the Raman lidar has a power-aperture product 500 times greater than that of the MPD.

1 Introduction

Water vapor (WV) is one of the fundamental thermodynamic variables that defines the state of the atmosphere and influences many important processes related to weather and climate. The importance of continuously monitoring lower tropospheric WV is underscored in the National Aeronautics and Space Administration (NASA) decadal survey (National Academies of Sciences, Engineering, and Medicine, 2018b) and in National Research Council (NRC, 2009, 2010, 2012) and National Academy of Sciences (National Academies of Sciences, Engineering, and Medicine, 2018a) reports. In particular, continuous range-resolved measurements of WV are needed at large scales to improve severe weather and precipitation predictions (Weckwerth et al., 1999; Wulfmeyer et al., 2015; Geerts et al., 2016; Jensen et al., 2016).

To fulfill this observational need, Montana State University (MSU) and the National Center for Atmospheric Research (NCAR) have developed a micropulse differential absorption lidar (MPD) that continuously measures range-resolved WV in the lower (150 m to 6 km) atmosphere (Spuler et al., 2015, 2021; NCAR/EOL MPD Team, 2020). The MPD is designed for unattended network deployment, using low-power, low-cost, high-reliability diode lasers that enable class 1M eye-safe transmitted lasers. The MPD employs the narrowband differential absorption lidar (DIAL) technique on a WV absorption line with low temperature sensitivity, whereby approximate knowledge of atmospheric temperature and pressure allow for a first-principles-based retrieval (Nehrir et al., 2009). An immediate benefit of the MPD is that it provides continuous WV measurements independent of radiosonde WV measurements, which can provide additional information to numerical weather prediction (NWP) data assimilation systems.

1.1 Problem statement – extending capability of MPD

The MPD has proven to produce accurate WV measurements up to 6 km altitude, depending on aerosol loading, clouds, and solar background. However, in high solar background conditions, MPD water vapor retrievals can be noisy as low as 2 km. The MPD's ability to produce precise measurements above 2 km and accurate higher-altitude WV measurements is impeded by the current standard WV DIAL retrieval methods. These methods are built upon a fundamental methodology that algebraically solves the WV variable through operations that exacerbate any random noise in the lidar observations (Marais et al., 2016). To suppress the noise the standard method applies low-pass filters on the algebraically computed WV using a Gaussian smoothing kernel or Savitzky–Golay filter (Schafer, 2011). The signal-to-noise ratio (SNR) of the WV measurements can be further improved upon by reducing the vertical and horizontal resolutions of the photon-counting observations. However, reducing the resolutions introduces systematic biases in the WV due to the nonlinearity of the single-scatter lidar equation. Smoothing operations are limited in their benefits because optimal averaging is generally localized to patches of correlated structure in a lidar profile. Invariably, parts of the profile are over-smoothed (resulting in biases), while other parts are under-smoothed (resulting in random error) (Hayman et al., 2020).

Figure 1 illustrates the shortcomings of low-pass filtering, where Fig. 1b and e show the standard method estimated WV image and profile obtained using low-pass filtering. The profile is juxtaposed with a coincident radiosonde (RS) WV profile. Although the bandwidth of the low-pass filter is relatively well suited for low altitude (≤ 4 km) WV measurements, Fig. 1i shows that the filter is not sufficient to meaningfully reduce the noise at higher altitudes. Furthermore, Fig. 1ii illustrate that the low-pass filter also over-smooths rapidly changing WV features such as embedded dry regions.

Figure 1 The (a) attenuated-backscatter cross section (atten. backs.) and the water vapor (WV) measurements of the (b) standard (STND) and (c) PTV-MPD methods. Panel (a) is not masked to show what atmospheric features have been masked out in (b) and (c) by the gray areas. The dashed vertical white lines in (b) and (c) show the launch time of the radiosonde (RS), where (e) and (f) show the RS WV profile. Panel (e) compares the standard method WV measurements with RS WV profile, whereas (f) shows the same comparison for the PTV-MPD method. The horizontal dashed lines in (e) to (f) show the altitude at which the RS was horizontally 5 km away from the MPD. The dashed boxes in (e) and (f) highlight regions that are discussed in the text.

1.2 Proposed approach to extending MPD capability

Advances have been made in denoising photon-counting medical images where the photon detection methodologies and forward modeling are similar to that in atmospheric lidar (Fessler, 2020; Oh et al., 2013; Harmany et al., 2012; Willett and Nowak, 2003; Ahn and Fessler, 2003). Basically, these methods quantitatively separate the image being estimated from the random noise by making a distinction between the (1) vertical and horizontal correlations among the pixels of the underlying image and (2) the uncorrelated photon-counting noise. The work of Xiao et al. (2020), Hayman and Spuler (2017), and Marais et al. (2016) have demonstrated that the medical denoising methods can be adopted and adapted to dramatically improve the inference of the extinction cross section from ground-based photon-counting high spectral resolution lidar (HSRL).

Inspired by the medical image denoising methods and the Poisson total variation (PTV) method for inferring extinction from HSRL observations (Xiao et al., 2020; Hayman and Spuler, 2017; Marais et al., 2016), we explore how the PTV method can be adapted for the MPD instrument to extend its capability for measuring WV. Figure 1c and f show an example WV measurement of the adapted PTV method, labeled as PTV-MPD, of the same scene and profile of Fig. 1c and e; using the range intervals indicated by Fig. 1i and ii, the PTV-MPD WV measurements are more accurate at higher altitudes and at the embedded dry region compared to the standard method.

The PTV-MPD method we present in this paper differs from PTV by three adaptations that are necessary for accurate WV measurements (Marais et al., 2016).

• Forward models for DIAL. The first adaptation of PTV requires that we develop MPD forward models to fit the estimated parameters on to the observed photon counts; while the objective is to estimate the WV, to accommodate the forward model, the attenuated backscatter is also an estimated product. To accurately model what the MPD observes, a convolution operator is included in the forward model that represents the oversampling of the laser pulse (Spuler et al., 2021); the MPD uses relatively long laser pulses to increase the SNR, with the trade-off of blurring the backscatter optical signal described by the standard single-scatter lidar equation (SSLE).

• Simultaneous inference. PTV for the HSRL infers the backscatter and extinction cross sections with a two-step approach in order to isolate geometric overlap calibration biases to the extinction cross section. However, with PTV-MPD a two-step approach is not necessary to infer the WV and attenuated backscatter since the DIAL employs a differential measurement technique, which decouples the WV measurement from the geometric overlap. Hence, the second adaptation of PTV is that PTV-MPD infers the WV and attenuated backscatter simultaneously instead of separately. The simultaneous inference approach produces WV estimates that are more faithful to the MPD observations compared to the two-step inference approach.

• Mitigate DIAL sensitivity to initial attenuated backscatter. The simultaneous inference of the WV and attenuated backscatter requires initial estimates of both the WV and attenuated backscatter. Our experiments show that PTV-MPD is sensitive to near-range inaccuracies in the initial attenuated-backscatter estimate, which induces inaccuracies in the near-range WV measurements. Therefore, the third adaptation involves making PTV-MPD more robust against inaccuracies in the initial attenuated backscatter.

1.3 Contributions

The immediate contributions of this paper are the

1. adaptation of the PTV method (i.e., PTV-MPD) for the DIAL technique and

2. the first rigorous validation of the PTV method using in situ measurements.

This work also serves as an example of how development of advanced signal processing can provide insights into improving hardware design and trades. By leveraging advances in signal processing, lidar hardware costs might be reduced whilst maintaining or improving the retrieval precision.

1.4 Outline, notation, and symbols

We introduce the MPD forward models in Sect. 2. Thereafter, in Sect. 3 we discuss the noise model of the MPD and how saturated photon counts are masked. The standard and PTV-MPD methods are discussed in Sects. 4 and 5. The paper ends with experiment results in Sect. 6 and the conclusion thereafter.

Table 1 Commonly used symbols.

Geophysical variables and the lidar forward models are written as matrices. When we introduce a geophysics-related variable we immediately indicate the units of the variable using the notation [⋅]; for example [W]. Table 1 lists the primary matrices and index variables used throughout this paper. The nth row and kth column element of a matrix is denoted, for example, by φ_n,k. We use the superscript index (ι) to indicate whether a matrix or vector is specific to lidar channel (i.e., wavelength).

Table 2 lists commonly used acronyms that are being used throughout this paper.

Table 2 Commonly used abbreviations.

2 The MPD forward models

To date NCAR has five experimental MPD instruments where each instrument transmits laser pulses at a rate of 7 to 8 kHz at two wavelengths sequentially. The first wavelength is tuned on a WV absorption line near 828.2 nm, whereas the second wavelength is off an WV absorption line near 828.3 nm (Spuler et al., 2015); these wavelengths are indexed by $\mathit{\iota }\in \mathit{\{}\mathrm{on},\mathrm{off}\mathit{\}}$. The MPDs transmit laser pulses that are longer than the sampling interval Δt=250 ns to increase the SNR of the WV measurements; the laser pulse duration was 1 µs during a Southern Great Plains (SGP) field campaign and has recently been decreased to 625 ns for the next generation of MPDs (Spuler et al., 2021). The integer ΔN≥0 quantifies how many vertical sampling intervals Δt, minus 1, fit in the duration of a laser pulse; for example, for the SGP data $\mathrm{\Delta }N=({\mathrm{10}}^{-\mathrm{6}}/\mathrm{\Delta }t)-\mathrm{1}=\mathrm{3}$.

The MPD instrument photon detector observes the weighted sum of the lidar equation and the laser pulse, since the laser pulse spans multiple range sampling intervals. We define the SSLE on the measurement range axis, which are denoted by ${\mathit{r}}_{n^{\prime }}$ where $n^{\prime }=\mathrm{1},\mathrm{2},\mathrm{\dots },N+\mathrm{\Delta }N$. The MPD samples the photon rates mid-duration of the laser pulse, and therefore the range axis of the observations is shifted relative to the measurement range axis. Specifically, the observation range axis is denoted by ${\stackrel{\mathrm{̃}}{\mathit{r}}}_n$, where $n=\mathrm{1},\mathrm{2},\mathrm{\dots },N$ indexes the observation range axis; the relation between the observation ${\stackrel{\mathrm{̃}}{\mathit{r}}}_n$ and measurement ${\mathit{r}}_{n^{\prime }}$ range axes is

$$\begin{array}{}\text{(1)}& {\stackrel{\mathrm{̃}}{\mathit{r}}}_n={\mathit{r}}_{n+\lceil \mathrm{\Delta }N/\mathrm{2}\rceil }.\end{array}$$

The SSLE is a function of the unknown WV ${\mathit{\phi }}_{n^{\prime },k}$ [g m⁻³] and uncalibrated attenuated backscatter ${\mathit{\chi }}_{n^{\prime },k}$ [m⁻¹ sr⁻¹], where the columns of the matrices are indexed by $k=\mathrm{1},\mathrm{2},\mathrm{\dots }K$. We write the SSLE as a matrix function since our attention is on inferring an image of the WV. The single-scatter lidar matrix function at measurement range ${\mathit{r}}_{n^{\prime }}$ for channel ι is denoted by

$$\begin{array}{}\text{(2)}& {\displaystyle }{\stackrel{\mathrm{̃}}{\mathbf{S}}}_{n^{\prime },k}^{\left(\mathit{\iota }\right)}(\mathit{\phi },\mathit{\chi })={\displaystyle \frac{{\mathbf{C}}_{n^{\prime },k}}{\mathrm{\Delta }t{\mathit{r}}_{n^{\prime }}^{\mathrm{2}}}}{\mathit{O}}_{n^{\prime }}{\mathit{\chi }}_{n^{\prime },k}\text{(3)}& {\displaystyle }\times \exp \left(-\mathrm{2}\mathrm{\Delta }r\sum _{n^{\prime \prime }=\mathrm{0}}^{n^{\prime }}{\mathit{\sigma }}_{n^{\prime \prime },k}^{\left(\mathit{\iota }\right)}{\mathit{\phi }}_{n^{\prime \prime },k}\right),\end{array}$$

and its unit is [W]. The calibration parameters are the range resolution $\mathrm{\Delta }r=c\mathrm{\Delta }t/\mathrm{2}\approx \mathrm{37.5}$ m, backscatter calibration constant ${\mathbf{C}}_{n^{\prime },k}$ [J sr m³], and geometric overlap ${\mathit{O}}_{n^{\prime }}$. The WV absorption [m² g⁻¹] of channel ι is denoted by ${\mathit{\sigma }}_{n^{\prime },k}^{\left(\mathit{\iota }\right)}$ and is pre-computed using (1) an assumed lapse rate with which the temperature and pressure profiles are approximated, (2) the mass of a water molecule per mole, and (3) the Avogadro number (Spuler et al., 2015; Nehrir et al., 2009).

The weighted sum of the SSLE with the laser pulse is modeled by

$$\begin{array}{}\text{(4)}& \mathrm{\Delta }t\sum _{n^{\prime }=n}^{n+\mathrm{\Delta }N}{\mathbf{A}}_{n,n^{\prime }}{\stackrel{\mathrm{̃}}{\mathbf{S}}}_{n^{\prime },k}^{\left(\mathit{\iota }\right)}(\mathit{\phi },\mathit{\chi }),\end{array}$$

where the matrix ${\mathbf{A}}_{n,n^{\prime }}$ models the laser energy distribution over the multiple range sampling intervals; the factor Δt models the integration of the observed photon rate by the photon detectors. The first ΔN+1 columns ($n^{\prime }=\mathrm{1},\mathrm{2},\mathrm{\dots },\mathrm{\Delta }N+\mathrm{1}$) of the first row (n=1) of the matrix ${\mathbf{A}}_{n,n^{\prime }}$ is the fractional laser energy distribution over the duration of the pulse such that the sum of the first row is equal to one, and the remaining N columns are equal to zero. The nth row of matrix ${\mathbf{A}}_{n,n^{\prime }}$ is the first row of ${\mathbf{A}}_{n,n^{\prime }}$ shifted circularly to the right n times. Hence, the MPD forward model at range ${\stackrel{\mathrm{̃}}{\mathit{r}}}_n$ for channel ι is defined by

$$\begin{array}{}\text{(5)}& {\mathbf{S}}_{n,k}^{\left(\mathit{\iota }\right)}\left(\mathit{\phi },\mathit{\chi }\right)={\mathit{U}}_k\mathrm{\Delta }t\left({\mathit{b}}_k^{\left(\mathit{\iota }\right)}+\sum _{n^{\prime }=n}^{n+\mathrm{\Delta }N}{\mathbf{A}}_{n,n^{\prime }}{\stackrel{\mathrm{̃}}{\mathbf{S}}}_{n^{\prime },k}^{\left(\mathit{\iota }\right)}(\mathit{\phi },\mathit{\chi })\right),\end{array}$$

where U_k is the number of laser shots per column index k. The vector ${\mathit{b}}_k^{\left(\mathit{\iota }\right)}$ is the dark and solar background photon rate [W] of channel ι.

The MPD instrument photon detectors are saturated over several range bins after each laser shot due to internal scattering in the co-axial optical configuration. Consequently, the first WV estimate starts at 125 m or 500 m depending on the hardware configuration. To model the unobserved WV from range 0 up to r₁, the lowest water absorption cross section ${\mathit{\sigma }}_{\mathrm{1},k}^{\left(\mathit{\iota }\right)}$ is the range axis sum of the absorption cross sections from range 0 up to r₁. Hence, the WV φ_1,k at range r₁ represents the range-average WV from range 0 up to r₁.

3 The MPD photon-counting noise model and masking

3.1 Photon-counting noise probability mass function

The photon-counting observations at range ${\stackrel{\mathrm{̃}}{\mathit{r}}}_n$ and profile k of the online and offline channels are denoted by ${\mathbf{Y}}_{n,k}^{\left(\mathit{\iota }\right)}$. Each photon count represents temporally accumulated counts of multiple laser shots as indicated by U_k. The noise of the photon counts are modeled by the Poisson probability mass function (PMF) if the corresponding photon rates are below the saturation limit of the MPD photon detectors (Donovan et al., 1993; Müller, 1973). The expected values of these unsaturated photon counts is modeled by the MPD forward model Eq. (5), and we assume that

$$\begin{array}{}\text{(6)}& \begin{array}{rl}{\mathbf{Y}}_{n,k}^{\left(\mathit{\iota }\right)}& \sim \mathrm{Poisson}\left({\mathbf{S}}_{n,k}^{\left(\mathit{\iota }\right)}\right)\\ & \equiv \exp \left(-{\mathbf{S}}_{n,k}^{\left(\mathit{\iota }\right)}+{\mathbf{Y}}_{n,k}^{\left(\mathit{\iota }\right)}{\log }_e{\mathbf{S}}_{n,k}^{\left(\mathit{\iota }\right)}-{\log }_e{\mathbf{Y}}_{n,k}^{\left(\mathit{\iota }\right)}\mathrm{!}\right)\end{array},\end{array}$$

where

$$\begin{array}{}\text{(7)}& {\mathbf{S}}_{n,k}^{\left(\mathit{\iota }\right)}\equiv {\mathbf{S}}_{n,k}^{\left(\mathit{\iota }\right)}\left(\mathit{\phi },\mathit{\chi }\right).\end{array}$$

3.2 Masking saturated photon counts

The instantaneous backscattered photon rates, corresponding to each laser shot, of clouds and precipitation can exceed the MPD photon detector saturation limit. Moreover, the accumulated photon counts ${\mathbf{Y}}_{n,k}^{\left(\mathit{\iota }\right)}$ consists of a combination of unsaturated and saturated photon counts within and at the edges of clouds and precipitation. Photon counts that are saturation contaminated cannot be accurately modeled by the Poisson PMF (Donovan et al., 1993). Hence, saturated photon counts and corresponding forward model pixels are excluded from our inference methodology using a mask matrix M_n,k.

The saturated photon count mask M_n,k has to be constructed via proxy data that indicate whether the instantaneous backscattered photon rates exceed the photon detector saturation limit. Since we know a priori that dense aerosol layers, clouds, and precipitation have large backscatter cross sections, the matrix M_n,k masks out these atmospheric features.

We used a sliding window standard deviation filter on the photon counts ${\mathbf{Y}}_{n,k}^{\left(\mathit{\iota }\right)}$, with a fixed threshold, to identify the large backscatter cross sections of dense aerosol layers, clouds, and precipitation; the threshold was set by qualitatively validating whether these large backscatter cross sections have been masked out.

3.3 Photon-counting noise negative log likelihood

The Poisson negative log likelihood (P-NLL) is used by the PTV-MPD method to quantify the fitting of the WV and attenuated backscatter onto the noisy observations via the MPD forward model Eq. (5). The P-NLL function of channel ι is denoted by

$$\begin{array}{}\text{(8)}& \begin{array}{rl}& {\mathcal{L}}_p\left({\mathbf{S}}^{\left(\mathit{\iota }\right)}\left(\mathit{\phi },\mathit{\chi }\right);{\mathbf{Y}}^{\left(\mathit{\iota }\right)}\right)=\sum _{n=\mathrm{1},k=\mathrm{1}}^{N,K}{\mathbf{M}}_{n,k}\\ & \times \left[p{\mathbf{S}}_{n,k}^{\left(\mathit{\iota }\right)}\left(\mathit{\phi },\mathit{\chi }\right)-{\mathbf{Y}}_{n,k}^{\left(\mathit{\iota }\right)}{\log }_e\left(p{\mathbf{S}}_{n,k}^{\left(\mathit{\iota }\right)}\left(\mathit{\phi },\mathit{\chi }\right)\right)\right],\end{array}\end{array}$$

where p is a normalization factor that is used in conjunction with the photon counts. The saturated photon counts are masked out with the matrix M_n,k. The factorial term of the Poisson PMF is omitted in Eq. (8) since it is a constant value.

4 The MPD standard method

A detailed discussion of the standard processing technique employed for MPD is provided in Spuler et al. (2021). Nonetheless, we will provide a brief overview of the method here as outlined in Alg. 1.

Algorithm 1 The MPD standard method.

Require:

The noisy observations ${\mathbf{Y}}_{n,k}^{\left(\mathrm{on}\right)}$ and ${\mathbf{Y}}_{n,k}^{\left(\mathrm{off}\right)}$.

1:

{Low pass filter background subtracted photon counts}

$$\begin{array}{}\text{(9)}& {\stackrel{\mathrm{̃}}{\mathbf{Y}}}_{n,k}^{\left(\mathit{\iota }\right)}=\mathrm{lowpass}-\mathrm{filter}\left({\mathbf{Y}}_{n,k}^{\left(\mathit{\iota }\right)}-{\mathit{U}}_k\mathrm{\Delta }t{\mathit{b}}_k^{\left(\mathit{\iota }\right)}\right)\end{array}$$

2:

{Compute differential optical depth}

$${\mathit{\tau }}_{n,k}=-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\log }_e\left({\displaystyle \frac{{\stackrel{\mathrm{̃}}{\mathbf{Y}}}_{n,k}^{\left(\mathrm{off}\right)}}{{\stackrel{\mathrm{̃}}{\mathbf{Y}}}_{n,k}^{\left(\mathrm{on}\right)}}}\right)$$

3:

{Compute the absolute water vapor}

$${\mathit{\tau }}_{\mathrm{0},k}=\mathbf{0},{\stackrel{\mathrm{̃}}{\mathit{\phi }}}_{n,k}={\displaystyle \frac{{\mathit{\tau }}_{n,k}-{\mathit{\tau }}_{n-\mathrm{1},k}}{{\mathit{\sigma }}_{n^{\prime \prime },k}^{\left(\mathrm{on}\right)}-{\mathit{\sigma }}_{n,k}^{\left(\mathrm{off}\right)}}}$$

4:

{Low pass filter the computed water vapor}

$${\widehat{\mathit{\phi }}}_{n,k}=\mathrm{lowpass}-\mathrm{filter}\left({\stackrel{\mathrm{̃}}{\mathit{\phi }}}_{n,k}\right)$$

5:

return  the water vapor estimate ${\widehat{\mathit{\phi }}}_{n,k}$.

The photon count profiles first undergo low-pass filtering to suppress the random noise (Hayman et al., 2020). The ratio of the online and offline channels cancels the attenuated backscatter in each forward model, which allows us to solve directly for the optical depth resulting from WV. At this stage, the WV estimate has residual noise that is low-pass-filtered with Gaussian kernels, with 5 to 10 min by 170 m bandwidths, to further reduce the noise in the observation. The smoothing kernel of the low-pass filter is fixed across the scene. Thus, the WV estimates are often over-smoothed in highly dynamic, high-SNR regions and under-smoothed at higher altitudes. Due to the direct nature of the inversion, there are also no constraints imposed on the retrieval, therefore nonphysical WV estimates are frequently obtained in noisy regions. The nonphysical WV values, in turn, can create problems in further downstream scientific analysis where nonphysical state variables present a challenge. Nonphysical quantities are generally not easily included in such analysis, but the selective omission or limiting of nonphysical noise will also create biases.

Finally, the lengths of the laser pulses are not accounted for in the standard method. Specifically, inclusion of the laser pulse in the forward model prevents the attenuated-backscatter term from directly canceling out when dividing the offline with the online channel, and a direct inversion becomes no longer possible. However, failing to account for the laser pulse length can create biases in some cases, for example at WV dry regions.

5 The PTV-MPD method

The PTV method, originally developed for photon-counting HSRL (Marais et al., 2016), is an adaptation of the SPIRAL method, which is a regularized maximum likelihood technique (Oh et al., 2013; Harmany et al., 2012); the technique and derivations of it have be applied in wide range of inverse problems in different domains, such as in medical imaging and astronomy (Roelofs et al., 2020; Harmany et al., 2012).

5.1 Overview of PTV

With PTV the assumptions are that

1. the photon-counting noise can be accurately modeled by the Poisson PMF in Eq. (6),

2. the expected value of the photon counts can be accurately modeled with a forward model (i.e., Eq. 5), and

3. the geophysical variable (i.e., WV) image that we want to estimate can be accurately approximated with a two-dimensional (2D) piecewise constant (PC) function.

An accurate noise model is important in inverse problems, since with Poisson noise the noise variance is signal dependent, and modeling the photon-counting noise as a Gaussian distribution will lead to inversion inaccuracies (Harmany et al., 2012). The 2D PC approximation induces spatial correlations on the estimated geophysical variable while preserving any discontinuities; such correlations are expected from the true unobserved geophysical variable of which WV and attenuated backscatter are examples.

By making a distinction between the random noise and the geophysical variable having spatial correlation, PTV is capable of separating the noise from the geophysical variable that is being estimated. PTV achieves this separation by (1) using the P-NLL from Eq. (8) to quantify how close of a fit the forward model is relative to the noisy observations and by (2) employing the total variation (TV) penalty function that regularizes the geophysical variable to be approximately 2D PC.

The approach used by PTV is formulated as a mathematical optimization framework where we search over all candidate geophysical variables and choose the geophysical variable that minimizes the sum of the (1) P-NLL composited with the forward model and (2) the TV penalty functions. The P-NLL composite is called the “loss function” and the loss summed with the penalty functions is called the “objective function”. The optimization framework is conceptually illustrated by the equation

$$\begin{array}{}\text{(10)}& \underset{\mathbf{X}}{min}\mathit{\left\{}\stackrel{\text{Objective function}}{\overbrace{\underset{\text{Loss function}}{\underbrace{\stackrel{\mathrm{̃}}{\mathrm{\ell }}(\mathbf{X};\mathbf{Y})}}+\underset{\text{Tuning parameter}}{\underbrace{\stackrel{\mathrm{̃}}{\mathit{\lambda }}}}\times \underset{\text{TV penalty function}}{\underbrace{\text{TV}\left(\mathbf{X}\right)}}}}\mathit{\right\}},\end{array}$$

where X ≡ geophysical variable, Y ≡ photon-counting observations, $\stackrel{\mathrm{̃}}{\mathrm{\ell }}(\mathbf{X};\mathbf{Y})$ ≡ P-NLL composited with forward model, e.g., Eq. (8).

The (anisotropic) TV penalty function is defined as (Harmany et al., 2012)

$$\begin{array}{}\text{(11)}& \text{TV}\left(\mathbf{X}\right)=\sum _{\begin{array}{c}n=\mathrm{1}\\ k=\mathrm{1}\end{array}}^{\begin{array}{c}N-\mathrm{1}\\ K\end{array}}\left|{\mathbf{X}}_{n,k}-{\mathbf{X}}_{n+\mathrm{1},k}\right|+\sum _{\begin{array}{c}n=\mathrm{1}\\ k=\mathrm{1}\end{array}}^{\begin{array}{c}N\\ K-\mathrm{1}\end{array}}\left|{\mathbf{X}}_{n,k}-{\mathbf{X}}_{n,k+\mathrm{1}}\right|,\end{array}$$

and is weighted with a tuning parameter $\stackrel{\mathrm{̃}}{\mathit{\lambda }}\ge \mathrm{0}$ that sets the degree to which the 2D PC regularization is promoted. The tuning parameter is quantitatively determined through a cross-validation methodology (Friedman et al., 2001, chap. 7)(Oh et al., 2013), and the parameter is not set by “expert opinion”.

5.2 PTV-MPD: adaptation required for MPD

We adopt the PTV method to infer the WV from the MPD observations, where we approximate both the WV and attenuated backscatter as 2D PC functions. In contrast, the standard method of Alg. 1 divides out the attenuated backscatter, which is not feasible with PTV-MPD since it is unclear how to model the noise of the ratio of background subtracted photon counts.

In order to apply PTV on MPD observations the following adaptations are required:

(A1)

use of the MPD forward models of Eq. (5),

(A2)

allowing for the simultaneous inference of the WV and attenuated backscatter.

(A3)

prevention of the degradation of the WV measurement due to inaccuracies in the initial attenuated-backscatter value.

Adaptations (A1) and (A2) are achieved by redefining the PTV loss and objective functions, as discussed in the next section. Specific to adaptation (A2) in the original HSRL formulation of PTV, (1) the backscatter cross section was computed from the denoised photon counts, and (2) the lidar ratio was inferred using the computed backscatter cross section. This two-step sequence isolates calibration biases due to the geometric overlap function from the inferred lidar ratio and the extinction cross section. With a DIAL instrument such as MPD, the geometric overlap does not induce a bias in the inferred WV, in part because both measurements use the same photon detector (Spuler et al., 2015, 2021). Therefore, the two-step sequence used for HSRL observations is not necessary.

Regarding adaptation (A3), the simultaneous inference of the WV and attenuated backscatter requires an initial value of the attenuated backscatter. During our investigation of adapting PTV for the MPD instrument we noticed that the adapted PTV method is sensitive to inaccuracies in the initial attenuated-backscatter value. Specifically, when computing the initial value of the attenuated backscatter algebraically from the photon-counting observations, the contributions to the inaccuracies in the initial attenuated backscatter originate from the small SNR at close ranges and the long laser pulses that convolve with the SSLE (see Eq. 4). Hence, for adaptation (A3) we propose that PTV-MPD employ a coarse-to-fine image resolution framework when inferring the WV (Marais and Willett, 2017). With the coarse-to-fine framework, we first infer the WV and attenuated backscatter at a coarse image resolution, which allows for reducing inaccuracies in the initial attenuated backscatter. In other words, at the coarse image resolution the long laser pulses are deconvolved from the initial attenuated backscatter and the SNR of close-range observations are implicitly increased. In a sequence from coarse-to-fine image resolution, we use the coarse-resolution WV and attenuated-backscatter estimates as initial values when inferring the finer image resolution WV and attenuated backscatter.

5.3 PTV-MPD method algorithmic details

Using the conceptual optimization framework in Eq. (10) as a reference, the PTV-MPD method is formulated as an optimization problem where the minimizers of the objective function are the estimates of the WV φ and attenuated backscatter χ. The objective function used by PTV-MPD is introduced in Sect. 5.3.1, which addresses adaptations (A1) and (A2); Sect. 5.3.1 serves as a preliminary section before elaborating on the coarse-to-fine image resolution framework. Minimizing an objective function requires initial values of the WV φ and attenuated backscatter χ and can have a significant impact on the estimate of these variables (Harmany et al., 2012; Kelley, 1999); to this effect, in Sect. 5.3.2 we discuss the initial values of the WV φ and attenuated backscatter χ. Finally, Sect. 5.3.3 discusses how the objective function is minimized through the coarse-to-fine image resolution framework that addresses adaption (A3).

The formulation of the PTV-MPD includes tuning parameters that are computed using a cross-validation (CV) methodology; the CV involves calculating a validation error that indicates the optimality of the tuning parameters and is discussed in Sect. 5.3.4. Related to the tuning parameters and validation error, we require a way to test the hypothesis that the coarse-to-fine image resolution framework will improve the accuracy of inferring the WV φ. To test the hypothesis, we compare the test errors from inferring the WV φ at only the finest image resolution versus from a coarse-to-fine image resolution. The test error computation is also discussed in Sect. 5.3.4.

Figure 2 Pictorial overview of Alg. 2 and its sub-algorithms (Algs. 3 and 4) for inferring the water vapor (WV) and uncalibrated attenuated-backscatter cross section (atten. backs.). The corresponding algorithm line numbers for each block in the flow diagrams are indicated by the vertical text Lx-y. The blue, violet, and red text is meant to highlight statistically independent datasets generated using Poisson thinning.

Figure 2 gives a broad pictorial overview of how each part of the PTV-MPD framework is implemented with modularized sub-algorithms; the purposes of the sub-algorithms are as follows.

Fig. 2a, Alg. 2

performs the necessary preparations to infer the WV φ and attenuated backscatter χ which includes computing the initial attenuated backscatter.

Fig. 2b, Alg. 3

employs a CV methodology to choose the optimum tuning parameters for inferring the WV φ and attenuated backscatter χ.

Fig. 2c, Alg. 4

infers the WV φ and attenuated backscatter χ for a fixed tuning parameter via a coarse-to-fine image resolution framework.

These algorithms are outlined in Sect. 5.3.5, which includes specific details about the SPIRAL method adaptations.

5.3.1 Preliminary for coarse-to-fine framework: the loss and objective functions

The PTV-MPD loss function differs from the PTV loss function in two respects; cf. Eq. (23) in Marais et al. (2016). First, the loss function of PTV-MPD is defined such that we infer the log _e of the attenuated backscatter denoted by ${\stackrel{\mathrm{̃}}{\mathit{\chi }}}_{n,k}\equiv {\log }_e{\mathit{\chi }}_{n,k}$, since previous work suggests that more accurate inference can be achieved by inferring the log _e of a linear variable with methods similar to PTV (Oh et al., 2013, 2014). Second, we define the PTV-MPD loss function as the sum of the normalized P-NLL function of each channel, where the normalization ensures that each P-NLL function is equally weighted; this weighting is done since the accumulated photon counts of the offline channels are generally larger than that of the online channel. Hence, the PTV-MPD loss function is defined by

$$\begin{array}{}\text{(12)}& \begin{array}{rl}& {\mathrm{\ell }}_p\left(\mathit{\phi },\stackrel{\mathrm{̃}}{\mathit{\chi }};{\mathbf{Y}}^{\left(\mathrm{on}\right)},{\mathbf{Y}}^{\left(\mathrm{off}\right)}\right)\\ & =\sum _{\mathit{\iota }\in \mathit{\{}\mathrm{on},\mathrm{off}\mathit{\}}}\underset{\text{channel }\mathit{\iota }\text{ P-NLL Eq. (8)}}{\underbrace{{\mathcal{L}}_p\left({\mathbf{S}}^{\left(\mathit{\iota }\right)}\left(\mathit{\phi },\exp \left(\stackrel{\mathrm{̃}}{\mathit{\chi }}\right)\right);{\mathbf{Y}}^{\left(\mathit{\iota }\right)}\right)}}\times \stackrel{\text{channel }\mathit{\iota }\text{ weight}}{\overbrace{{\mathit{\omega }}^{\left(\mathit{\iota }\right)}}}\end{array},\end{array}$$

which is the weighted sum of the P-NLL of each channel where each P-NLL is normalized by the root sum square (i.e., Frobenius norm) of the photon counts that have not been masked out. The Frobenius norm of the mask photon counts is defined by

$$\begin{array}{}\text{(13)}& {\mathit{\omega }}^{\left(\mathit{\iota }\right)}={∥{\mathbf{MY}}^{\left(\mathit{\iota }\right)}∥}_F^{-\mathrm{1}}={\left[\sum _{\begin{array}{c}n=\mathrm{1},\\ k=\mathrm{1}\end{array}}^{N,K}{\left({\mathbf{M}}_{n,k}{\mathbf{Y}}_{n,k}^{\left(\mathit{\iota }\right)}\right)}^{\mathrm{2}}\right]}^{-\mathrm{1}/\mathrm{2}}.\end{array}$$

The objective function that is minimized by PTV-MPD is defined as

$$\begin{array}{}\text{(14)}& \begin{array}{rl}& {\stackrel{\mathrm{̃}}{F}}_{p,{\mathit{\lambda }}_{\mathrm{w}},{\mathit{\lambda }}_{\mathrm{a}}}\left(\mathit{\phi },\stackrel{\mathrm{̃}}{\mathit{\chi }};{\mathbf{Y}}^{\left(\mathrm{on}\right)},{\mathbf{Y}}^{\left(\mathrm{off}\right)}\right)\\ & \equiv \underset{\text{Loss function, Eq. (12)}}{\underbrace{{\mathrm{\ell }}_p\left(\mathit{\phi },\stackrel{\mathrm{̃}}{\mathit{\chi }};{\mathbf{Y}}^{\left(\mathrm{on}\right)},{\mathbf{Y}}^{\left(\mathrm{off}\right)}\right)}}+\stackrel{\text{Regularize WV}}{\overbrace{{\mathit{\lambda }}_{\mathrm{w}}\text{TV}\left(\mathit{\phi }\right)}}\\ & +\underset{\text{Regularize atten. backs.}}{\underbrace{{\mathit{\lambda }}_{\mathrm{a}}\text{TV}\left(\stackrel{\mathrm{̃}}{\mathit{\chi }}\right),}},\end{array}\end{array}$$

which is parameterized with the WV φ and attenuated backscatter $\stackrel{\mathrm{̃}}{\mathit{\chi }}$ and their respective tuning parameters λ_w and λ_a. The WV λ_w and attenuated backscatter λ_a tuning parameters set the degree to which these parameters should be regularized to be 2D PC. While minimizing the objective function (14), the WV φ is constrained to be non-negative, which is denoted by the set 𝒲. The objective function (14) is minimized using the alternating minimization method in conjunction with adaptations of SPIRAL (Beck and Tetruashvili, 2013; Harmany et al., 2012); the alternating minimization is elaborated in Sect. 5.3.5 and outlined by Alg. 4.

5.3.2 Initial values of objective function

In Appendix B we show that for low-SNR regions of the photon-counting images there can be multiple estimates of the WV φ or attenuated backscatter $\stackrel{\mathrm{̃}}{\mathit{\chi }}$ that minimizes the objective function (14). The dark and solar background photon rates b^(ι) are the primary factors that determine the domains over which the objective function has a unique minimizer. Consequently, we expect that PTV-MPD will require accurate initial values of the WV φ and attenuated backscatter $\stackrel{\mathrm{̃}}{\mathit{\chi }}$ whenever the photon count SNRs are dominated by the background rate b^(ι) such as at close range and in high-altitude regions.

It will be advantageous for the MPD to make WV measurements that are statistically independent of RS for the purpose of providing additional statistically independent information for NWP data assimilation. Thus, in this paper we set the initial WV to zero [g m⁻³]; this WV initial value reflects that we assume no a priori information about the WV other than (1) the assumed lapse rate¹ of the atmosphere and (2) the assumption that the WV can be accurately approximated with 2D PC functions. This initial WV value is denoted by ${\widehat{\mathit{\phi }}}^{\text{init}}$.

The offline channel is less sensitive to WV absorption compared to the online channel. Hence, the initial value of the attenuated backscatter, denoted by ${\widehat{\stackrel{\mathrm{̃}}{\mathit{\chi }}}}^{\text{init}}$, is computed from the offline channel with the initial WV ${\widehat{\mathit{\phi }}}^{\text{init}}$ via

$$\begin{array}{}\text{(15)}& {\widehat{\stackrel{\mathrm{̃}}{\mathit{\chi }}}}_{n^{\prime },k}^{\text{init}}={\log }_e\left({\displaystyle \frac{{\sum }_{n=\mathrm{1}}^N{\mathbf{A}}_{n^{\prime },n}^{\mathrm{T}}\left({\mathbf{Y}}_{n,k}^{\left(\mathrm{off}\right)}-{\mathit{U}}_k\mathrm{\Delta }t{\mathit{b}}_k^{\left(\mathit{\iota }\right)}\right)}{{\mathit{U}}_k\mathrm{\Delta }t{\stackrel{\mathrm{̃}}{\mathbf{S}}}_{n^{\prime },k}^{\left(\mathrm{off}\right)}({\widehat{\mathit{\phi }}}^{\text{init}},\mathbf{1})}}\right).\end{array}$$

The matrix A^T maps the observations to the measurement space; the matrix is the adjoint (i.e., transpose) of the laser energy distribution matrix A. The vector U is the number of accumulated laser shots, and Δt models the integration of the observed photon rate by the photon detectors (see Sect. 2).

5.3.3 The coarse-to-fine inference framework

With the coarse-to-fine image resolution inference framework, we start with estimating a coarse-resolution WV image and use the coarse-resolution estimate as an initial WV value at a finer image resolution; this technique has proven useful to more accurately denoising and inverting images for low-SNR photon-limited application (Marais and Willett, 2017; Azzari and Foi, 2017). The coarse image resolution estimation is akin to downsampling the noisy observations to increase the SNR, though the difference is that with the proposed methodology the downsampling in the coarse-to-fine framework takes into account the MPD forward model (5).

Following the approaches delineated in Marais and Willett (2017) and Azzari and Foi (2017), we define the following variables and linear operators.

1. A series of values $\stackrel{\mathrm{‾}}{h}=h_{\mathrm{1}}>h_{\mathrm{2}}>\mathrm{\dots }>h_L=\mathrm{1}$ that represents the downsampling factors of the coarse-to-fine image resolutions at which the WV will be inferred; $\stackrel{\mathrm{‾}}{h}$ represents the coarsest image resolution.

2. The downsampling operator is denoted by ${\mathcal{D}}_h^{\downarrow }$ and its corresponding upsampling operator ${\mathcal{D}}_h^{\uparrow }$.

The downsampled image ${\mathcal{D}}_h^{\downarrow }\mathit{\phi }$ will have approximately $N/h$ rows and $K/h$ columns, depending on how the downsample operator handles boundary pixels. The coarse-resolution WV image is denoted by ϕ_h, and its relationship with the finest image resolution WV image φ is governed by the upsampling operator ${\mathcal{D}}_h^{\uparrow }$ where $\mathit{\phi }={\mathcal{D}}_h^{\uparrow }{\mathit{\varphi }}_h$. The corresponding attenuated backscatter is denoted by ψ_h and is not downsampled like the WV; the subscript h in ψ_h indicates that it is paired with the downsampled WV ϕ_h estimate.

Let $l=\mathrm{1},\mathrm{2},\mathrm{\dots },L$ represent the coarse-to-fine image resolution index. The coarse-to-fine process repeats the following steps for each l until the finest image resolution WV estimate is obtained.

Step 1. Define the current image resolution WV estimate ${\mathit{\varphi }}_{h_l}$ and its paired attenuated-backscatter estimate ${\mathit{\psi }}_{h_l}$. In addition, define the current image resolution objective function.

$$\begin{array}{}\text{(16)}& \begin{array}{rl}& F_{p,{\mathit{\lambda }}_{\mathrm{w}},{\mathit{\lambda }}_{\mathrm{a}}}\left({\mathcal{D}}_{h_l}^{\uparrow }{\mathit{\varphi }}_{h_l},{\mathit{\psi }}_{h_l};{\mathbf{Y}}^{\left(\mathrm{on}\right)},{\mathbf{Y}}^{\left(\mathrm{off}\right)}\right)\\ & \equiv {\mathrm{\ell }}_p\left({\mathcal{D}}_{h_l}^{\uparrow }{\mathit{\varphi }}_{h_l},{\mathit{\psi }}_{h_l};{\mathbf{Y}}^{\left(\mathrm{on}\right)},{\mathbf{Y}}^{\left(\mathrm{off}\right)}\right)+{\mathit{\lambda }}_{\mathrm{w}}\text{TV}\left({\mathit{\varphi }}_{h_l}\right)\\ & +{\mathit{\lambda }}_{\mathrm{a}}\text{TV}\left({\mathit{\psi }}_{h_l}\right).\end{array}\end{array}$$

Step 2. If l=1, use the provided initial value of the WV ${\mathit{\varphi }}_{h_l}$ and attenuated backscatter ${\mathit{\psi }}_{h_l}$ when minimizing the objective function (16). Otherwise, use the previous image resolution estimates as initial values for the current image resolution. Specifically, the initial value for the WV is ${\mathcal{D}}_{h_l}^{\downarrow }{\mathcal{D}}_{h_{l-\mathrm{1}}}^{\uparrow }{\widehat{\mathit{\varphi }}}_{h_{l-\mathrm{1}}}$, and for the attenuated backscatter it is ${\widehat{\mathit{\psi }}}_{h_{l-\mathrm{1}}}$.

Step 3. Estimate the current image resolution WV ${\mathit{\varphi }}_{h_l}$ and its paired attenuated backscatter ${\mathit{\psi }}_{h_l}$ by minimizing the objective function (16):

$$\begin{array}{}\text{(17)}& \begin{array}{rl}& {\widehat{\mathit{\varphi }}}_{h_l},{\widehat{\mathit{\psi }}}_{h_l}\\ & =\arg \underset{{\mathit{\varphi }}_{h_l}\in \mathcal{W},{\mathit{\psi }}_{h_l}}{min}{F}_{p,{\mathit{\lambda }}_{\mathrm{w}},{\mathit{\lambda }}_{\mathrm{a}}}\\ & \left(\stackrel{\text{Fitting is performed on native resolution...}}{\overbrace{\underset{\text{...by upsampling coarse resolution WV image.}}{\underbrace{{\mathcal{D}}_{h_l}^{\uparrow }}}{\mathit{\varphi }}_{h_l}}},{\mathit{\psi }}_{h_l};{\mathbf{Y}}^{\left(\mathrm{on}\right)},{\mathbf{Y}}^{\left(\mathrm{off}\right)}\right)\end{array},\end{array}$$

where 𝒲 is the set of non-negative numbers.

Final step. If l=L, the final WV and attenuated-backscatter estimates are denoted by ${\widehat{\mathit{\phi }}}_{\stackrel{\mathrm{‾}}{h},{\mathit{\lambda }}_{\mathrm{w}},{\mathit{\lambda }}_{\mathrm{a}}}$ and ${\widehat{\stackrel{\mathrm{̃}}{\mathit{\chi }}}}_{\stackrel{\mathrm{‾}}{h},{\mathit{\lambda }}_{\mathrm{w}},{\mathit{\lambda }}_{\mathrm{a}}}$ meaning

$${\widehat{\mathit{\phi }}}_{\stackrel{\mathrm{‾}}{h},{\mathit{\lambda }}_{\mathrm{w}},{\mathit{\lambda }}_{\mathrm{a}}}={\widehat{\mathit{\varphi }}}_{h_L}$$

and

$$\begin{array}{}\text{(18)}& {\widehat{\stackrel{\mathrm{̃}}{\mathit{\chi }}}}_{\stackrel{\mathrm{‾}}{h},{\mathit{\lambda }}_{\mathrm{w}},{\mathit{\lambda }}_{\mathrm{a}}}={\widehat{\mathit{\psi }}}_{h_L}.\end{array}$$

The subscripts $\stackrel{\mathrm{‾}}{h}$, λ_w and λ_a indicate that the estimates of the WV φ and attenuated backscatter $\stackrel{\mathrm{̃}}{\mathit{\chi }}$ are the specific to coarse-to-fine image resolution configuration parameterized by $\stackrel{\mathrm{‾}}{h}$ and the tuning parameters λ_w and λ_a.

In practice, the minimization of the objective function (16) is done by alternating minimization between the WV and attenuated backscatter as indicated in Fig. 2c and outlined by Alg. 4.

5.3.4 Computing validation error to choose optimum tuning parameters and computing test error

The validation and test errors are computed by basically comparing the WV φ and attenuated backscatter $\stackrel{\mathrm{̃}}{\mathit{\chi }}$ estimates against noisy observations that are statistically independent of these estimates. This is achieved by thinning the non-masked offline Y^(off) and online Y^(on) photon counts through the Poisson thinning technique (Marais et al., 2016; Oh et al., 2013; Hayman et al., 2020); these thinned photon counts are statistically independent of each other and are Poisson distributed (Cinlar, 2013, chap. 4). Specifically, for each non-masked pixel in a photon-counting image Y^(ι), the Poisson thinning technique randomly samples the individual photon counts, and the sampled photon counts are placed in three photon-counting images. We call these three images the training ${\mathbf{Y}}_{\text{trn}}^{\left(\mathit{\iota }\right)}$, validation ${\mathbf{Y}}_{\text{vld}}^{\left(\mathit{\iota }\right)}$, and test ${\mathbf{Y}}_{\text{tst}}^{\left(\mathit{\iota }\right)}$ photon-counting images. The expected value of the training, validation, and test photon counts relative to the original photon counts are expressed by the scalars $p_{\text{trn}}+p_{\text{vld}}+p_{\text{tst}}=\mathrm{1}$, such that $\mathbb{E}\left[{\mathbf{Y}}_{\text{trn}}^{\left(\mathit{\iota }\right)}\right]=p_{\text{trn}}\mathbb{E}\left[{\mathbf{Y}}^{\left(\mathit{\iota }\right)}\right]$, where $\mathbb{E}\left[{\mathbf{Y}}_{\text{vld}}^{\left(\mathit{\iota }\right)}\right]$ and $\mathbb{E}\left[{\mathbf{Y}}_{\text{tst}}^{\left(\mathit{\iota }\right)}\right]$ are similarly defined.

Training photon counts

The training photon counts ${\mathbf{Y}}_{\text{trn}}^{\left(\mathit{\iota }\right)}$ are used to infer the WV ${\widehat{\mathit{\phi }}}_{\stackrel{\mathrm{‾}}{h},{\mathit{\lambda }}_{\mathrm{w}},{\mathit{\lambda }}_{\mathrm{a}}}$ and attenuated backscatter ${\widehat{\stackrel{\mathrm{̃}}{\mathit{\chi }}}}_{\stackrel{\mathrm{‾}}{h},{\mathit{\lambda }}_{\mathrm{w}},{\mathit{\lambda }}_{\mathrm{a}}}$ for specific WV λ_w and attenuated backscatter λ_a tuning parameters and coarse-to-fine configuration $\stackrel{\mathrm{‾}}{h}$. Specifically, when using the training photon counts the objective function in Eq. (16) will be equal to

$$\begin{array}{}\text{(19)}& F_{p_{\text{trn}},{\mathit{\lambda }}_{\mathrm{w}},{\mathit{\lambda }}_{\mathrm{a}}}\left({\mathcal{D}}_{h_l}^{\uparrow }{\mathit{\varphi }}_{h_l},{\mathit{\psi }}_{h_l};{\mathbf{Y}}_{\text{trn}}^{\left(\mathrm{on}\right)},{\mathbf{Y}}_{\text{trn}}^{\left(\mathrm{off}\right)}\right),\end{array}$$

where p in Eq. (16) has been replaced with p_trn.

Validation photon counts

The validation error of the estimates per tuning parameter is computed with the P-NLL (8) and the forward model (5) and is denoted by

$$\begin{array}{}\text{(20)}& \begin{array}{rl}& {\mathrm{err}}_{\mathrm{vld},{\mathit{\lambda }}_{\mathrm{w}},{\mathit{\lambda }}_{\mathrm{a}}}\\ & ={\mathcal{L}}_{p_{\text{vld}}}\left({\mathbf{S}}^{\left(\mathrm{on}\right)}\left({\widehat{\mathit{\phi }}}_{\stackrel{\mathrm{‾}}{h},{\mathit{\lambda }}_{\mathrm{w}},{\mathit{\lambda }}_{\mathrm{a}}},\exp \left({\widehat{\stackrel{\mathrm{̃}}{\mathit{\chi }}}}_{\stackrel{\mathrm{‾}}{h},{\mathit{\lambda }}_{\mathrm{w}},{\mathit{\lambda }}_{\mathrm{a}}}\right)\right);{\mathbf{Y}}_{\text{vld}}^{\left(\mathrm{on}\right)}\right)\\ & +{\mathcal{L}}_{p_{\text{vld}}}\left({\mathbf{S}}^{\left(\mathrm{off}\right)}\left({\widehat{\mathit{\phi }}}_{\stackrel{\mathrm{‾}}{h},{\mathit{\lambda }}_{\mathrm{w}},{\mathit{\lambda }}_{\mathrm{a}}},\exp \left({\widehat{\stackrel{\mathrm{̃}}{\mathit{\chi }}}}_{\stackrel{\mathrm{‾}}{h},{\mathit{\lambda }}_{\mathrm{w}},{\mathit{\lambda }}_{\mathrm{a}}}\right)\right);{\mathbf{Y}}_{\text{vld}}^{\left(\mathrm{off}\right)}\right).\end{array}\end{array}$$

The optimum tuning parameters corresponds to the smallest validation error, and this give us the WV and attenuated backscatter estimates

$$\begin{array}{}\text{(21)}& {\widehat{\mathit{\phi }}}_{\stackrel{\mathrm{‾}}{h}},{\widehat{\stackrel{\mathrm{̃}}{\mathit{\chi }}}}_{\stackrel{\mathrm{‾}}{h}}=\arg \underset{{\mathit{\lambda }}_{\mathrm{w}},{\mathit{\lambda }}_{\mathrm{a}}}{min}{\mathrm{err}}_{\mathrm{vld},{\mathit{\lambda }}_{\mathrm{w}},{\mathit{\lambda }}_{\mathrm{a}}}.\end{array}$$

Test photon counts

The test error for the estimates ${\widehat{\mathit{\phi }}}_{\stackrel{\mathrm{‾}}{h}}$ and ${\widehat{\stackrel{\mathrm{̃}}{\mathit{\chi }}}}_{\stackrel{\mathrm{‾}}{h}}$ are also computed with the P-NLL (8) and the forward model (5) and is denoted by

$$\begin{array}{}\text{(22)}& \begin{array}{rl}{\mathrm{err}}_{\mathrm{tst},\stackrel{\mathrm{‾}}{h}}& ={\mathcal{L}}_{p_{\text{tst}}}\left({\mathbf{S}}^{\left(\mathrm{on}\right)}\left({\widehat{\mathit{\phi }}}_{\stackrel{\mathrm{‾}}{h}},\exp \left({\widehat{\stackrel{\mathrm{̃}}{\mathit{\chi }}}}_{\stackrel{\mathrm{‾}}{h}}\right)\right);{\mathbf{Y}}_{\text{tst}}^{\left(\mathrm{on}\right)}\right)\\ & +{\mathcal{L}}_{p_{\text{tst}}}\left({\mathbf{S}}^{\left(\mathrm{off}\right)}\left({\widehat{\mathit{\phi }}}_{\stackrel{\mathrm{‾}}{h}},\exp \left({\widehat{\stackrel{\mathrm{̃}}{\mathit{\chi }}}}_{\stackrel{\mathrm{‾}}{h}}\right)\right);{\mathbf{Y}}_{\text{tst}}^{\left(\mathrm{off}\right)}\right).\end{array}\end{array}$$

To test the hypotheses as to whether the coarse-to-fine image resolution framework does improve the WV measurements, we evaluate the inequality

$$\begin{array}{}\text{(23)}& {\mathrm{err}}_{\mathrm{tst},\stackrel{\mathrm{‾}}{h}=\mathrm{9}}<{\mathrm{err}}_{\mathrm{tst},\stackrel{\mathrm{‾}}{h}=\mathrm{1}}.\end{array}$$

5.3.5 The PTV-MPD algorithm

Algorithm 2​​​​​​​ and its sub-algorithms (Algs. 3 and 4) outline the PTV-MPD algorithms as described by Fig. 2. Algorithm 2 takes as input the photon counts and the coarse image resolution at which PTV-MPD will start inferring the WV; the starting image resolution is expressed as a downsampling factor $\stackrel{\mathrm{‾}}{h}\ge \mathrm{1}$. Algorithm 3 outlines the CV methodology, and Alg. 4 outlines the coarse-to-fine WV inference image resolution framework.

The alternating minimization of the objective function (17) (Beck and Tetruashvili, 2013), employed by Alg. 4, is achieved by adaptations of SPIRAL where each adaptation is specific to the WV and attenuated backscatter variable (Oh et al., 2013; Harmany et al., 2012). The adaptation of SPIRAL includes replacing the original objective function and gradient matrix of the loss function; Appendix A presents the loss function gradient matrix along with the Jacobian matrix of the MPD forward model. The stopping criteria that is employed to stop the alternating minimization iterations in Alg. 4 uses the relative distance between consecutive estimates at the specific coarse image resolution. In more detail, if between iterations t+1 and t the relative distance, i.e.,

$$\begin{array}{}\text{(24)}& {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\displaystyle \frac{{∥{\mathit{\varphi }}^{(t+\mathrm{1})}-{\mathit{\varphi }}^{\left(t\right)}∥}_F}{{∥{\mathit{\varphi }}^{(t+\mathrm{1})}∥}_F}}+{\displaystyle \frac{{∥{\mathit{\psi }}^{(t+\mathrm{1})}-{\mathit{\psi }}^{\left(t\right)}∥}_F}{{∥{\mathit{\psi }}^{(t+\mathrm{1})}∥}_F}}\right),\end{array}$$

is less than 10⁻⁵ then the loop in Alg. 4 terminates, where

$$\begin{array}{}\text{(25)}& {∥\mathit{\varphi }∥}_F={\left(\sum _{n=\mathrm{1},k=\mathrm{1}}^{N,K}{\mathit{\varphi }}_{n,k}^{\mathrm{2}}\right)}^{\mathrm{1}/\mathrm{2}}.\end{array}$$

Algorithm 2 The PTV-MPD method.

Require:

(1) The noisy observations Y^(on) and Y^(off), (2) (mandatory) the coarsest resolution downsampling factor $\stackrel{\mathrm{‾}}{h}$, (3) (optional) initial values of the WV ${\widehat{\mathit{\phi }}}^{\text{init}}$, (4) (optional) tuning parameters for WV λ_w and atten. backs. λ_a.

1:

if λ_a and λ_w are not provided then

2:

${\mathrm{\Lambda }}_{\stackrel{\mathrm{̃}}{\mathit{\chi }}}={\mathrm{\Lambda }}_{\stackrel{\mathrm{̃}}{\mathit{\phi }}}={\left\{{\mathrm{10}}^{-\mathrm{2}+i/\mathrm{11}}\right\}}_{i=\mathrm{0}}^{\mathrm{11}}=\left\{{\mathrm{10}}^{-\mathrm{2}},{\mathrm{10}}^{-\mathrm{1.91}},\mathrm{\dots },{\mathrm{10}}^{\mathrm{2}}\right\}$

3:

else

4:

${\mathrm{\Lambda }}_{\stackrel{\mathrm{̃}}{\mathit{\chi }}}=\left\{{\mathit{\lambda }}_{\mathrm{a}}\right\}$ and ${\mathrm{\Lambda }}_{\stackrel{\mathrm{̃}}{\mathit{\phi }}}=\left\{{\mathit{\lambda }}_{\mathrm{w}}\right\}$.

5:

end if

6:

If WV ${\widehat{\mathit{\phi }}}^{\text{init}}$ is not provided, set ${\widehat{\mathit{\phi }}}^{\text{init}}=\mathbf{0}$

7:

{Split photon counts in training, validation and test counts}

8:

${\mathbf{Y}}_{\text{trn}}^{\left(\mathit{\iota }\right)},{\mathbf{Y}}_{\text{vld}}^{\left(\mathit{\iota }\right)},{\mathbf{Y}}_{\text{tst}}^{\left(\mathit{\iota }\right)}=$ Poisson Thinning (Y^(ι))

9:

{Compute the atten. backs. initial value}

10:

${\widehat{\stackrel{\mathrm{̃}}{\mathit{\chi }}}}^{\text{init}}$ = plug initial WV ${\widehat{\mathit{\phi }}}^{\text{init}}$ into Eq. (15) using scaled training photon counts ${\mathbf{Y}}_{\text{trn}}^{\left(\mathrm{off}\right)}/p_{\text{trn}}$

11:

{Infer the WV and atten. backs.}

12:

$$\begin{array}{}\text{(26)}& \begin{array}{rl}{\widehat{\mathit{\phi }}}_{\stackrel{\mathrm{‾}}{h}},{\widehat{\stackrel{\mathrm{̃}}{\mathit{\chi }}}}_{\stackrel{\mathrm{‾}}{h}}=\text{Algorithm 3}(& {\widehat{\mathit{\phi }}}^{\text{init}},{\widehat{\stackrel{\mathrm{̃}}{\mathit{\chi }}}}^{\text{init}},\stackrel{\mathrm{‾}}{h},{\mathrm{\Lambda }}_{\stackrel{\mathrm{̃}}{\mathit{\chi }}},{\mathrm{\Lambda }}_{\stackrel{\mathrm{̃}}{\mathit{\phi }}},\\ & {\mathbf{Y}}_{\text{trn}}^{\left(\mathrm{on}\right)},{\mathbf{Y}}_{\text{vld}}^{\left(\mathrm{on}\right)},{\mathbf{Y}}_{\text{trn}}^{\left(\mathrm{off}\right)},{\mathbf{Y}}_{\text{vld}}^{\left(\mathrm{off}\right)})\end{array}\end{array}$$

13:

{Compute the atten. backs.}

14:

${\widehat{\mathit{\chi }}}_{\stackrel{\mathrm{‾}}{h}}=\exp \left({\widehat{\stackrel{\mathrm{̃}}{\mathit{\chi }}}}_{\stackrel{\mathrm{‾}}{h}}\right)$

15:

Compute test error ${\mathrm{err}}_{\mathrm{tst},\stackrel{\mathrm{‾}}{h}}$ via Eq. (22)

16:

return  ${\widehat{\mathit{\phi }}}_{\stackrel{\mathrm{‾}}{h}},{\widehat{\mathit{\chi }}}_{\stackrel{\mathrm{‾}}{h}},{\mathrm{err}}_{\mathrm{tst},\stackrel{\mathrm{‾}}{h}}$

Algorithm 3 Cross-validation methodology in choosing optimum tuning parameters when inferring WV φ and attenuated backscatter χ.

Require:

(1) WV ${\widehat{\mathit{\phi }}}^{\text{init}}$ and atten. backs. ${\widehat{\stackrel{\mathrm{̃}}{\mathit{\chi }}}}^{\text{init}}$ initial values, (2) the coarsest resolution downsampling factor $\stackrel{\mathrm{‾}}{h}$, (3) tuning parameter sets ${\mathrm{\Lambda }}_{\stackrel{\mathrm{̃}}{\mathit{\chi }}}$ and ${\mathrm{\Lambda }}_{\stackrel{\mathrm{̃}}{\mathit{\phi }}}$, (4) training ${\mathbf{Y}}_{\text{trn}}^{\left(\mathrm{on}\right)},{\mathbf{Y}}_{\text{trn}}^{\left(\mathrm{off}\right)}$ and validation ${\mathbf{Y}}_{\text{vld}}^{\left(\mathrm{on}\right)},{\mathbf{Y}}_{\text{vld}}^{\left(\mathrm{off}\right)}$ photon counts.

1:

{Choose tuning parameters via cross-validation}

2:

for $({\mathit{\lambda }}_{\mathrm{w}},{\mathit{\lambda }}_{\mathrm{a}})\in {\mathrm{\Lambda }}_{\stackrel{\mathrm{̃}}{\mathit{\chi }}}\times {\mathrm{\Lambda }}_{\stackrel{\mathrm{̃}}{\mathit{\phi }}}$ do

3:

{Infer WV and atten. backs. using redefined objective function}

4:

${\widehat{\mathit{\phi }}}_{\stackrel{\mathrm{‾}}{h},{\mathit{\lambda }}_{\mathrm{w}},{\mathit{\lambda }}_{\mathrm{a}}},{\widehat{\stackrel{\mathrm{̃}}{\mathit{\chi }}}}_{\stackrel{\mathrm{‾}}{h},{\mathit{\lambda }}_{\mathrm{w}},{\mathit{\lambda }}_{\mathrm{a}}}$ =

Algorithm 4 $\left({\widehat{\mathit{\phi }}}^{\text{init}}{\widehat{\stackrel{\mathrm{̃}}{\mathit{\chi }}}}^{\text{init}},{\mathit{\lambda }}_{\mathrm{w}},{\mathit{\lambda }}_{\mathrm{a}},{\mathbf{Y}}_{\text{trn}}^{\left(\mathrm{on}\right)}{\mathbf{Y}}_{\text{trn}}^{\left(\mathrm{off}\right)},\stackrel{\mathrm{‾}}{h}\right)$

5:

Compute the validation error ${\mathrm{err}}_{\mathrm{vld},{\mathit{\lambda }}_{\mathrm{w}},{\mathit{\lambda }}_{\mathrm{a}}}$ via Eq. (20)

6:

end for

7:

{Choose tuning parameter that has smallest validation error}

$$\begin{array}{}\text{(27)}& {\mathit{\lambda }}_{\mathrm{w}}^{*},{\mathit{\lambda }}_{\mathrm{a}}^{*}=\arg \underset{{\mathit{\lambda }}_{\mathrm{w}},{\mathit{\lambda }}_{\mathrm{a}}}{min}\left\{{\mathrm{err}}_{\mathrm{vld},{\mathit{\lambda }}_{\mathrm{w}},{\mathit{\lambda }}_{\mathrm{a}}}\right\}\end{array}$$

8:

{With chosen tuning parameters select estimates}

9:

${\widehat{\mathit{\phi }}}_{\stackrel{\mathrm{‾}}{h}}={\widehat{\mathit{\phi }}}_{\stackrel{\mathrm{‾}}{h},{\mathit{\lambda }}_{\mathrm{w}}^{*},{\mathit{\lambda }}_{\mathrm{a}}^{*}}$, ${\widehat{\stackrel{\mathrm{̃}}{\mathit{\chi }}}}_{\stackrel{\mathrm{‾}}{h}}={\widehat{\stackrel{\mathrm{̃}}{\mathit{\chi }}}}_{\stackrel{\mathrm{‾}}{h},{\mathit{\lambda }}_{\mathrm{w}}^{*},{\mathit{\lambda }}_{\mathrm{a}}^{*}}$

10:

return  ${\widehat{\mathit{\phi }}}_{\stackrel{\mathrm{‾}}{h}},{\widehat{\stackrel{\mathrm{̃}}{\mathit{\chi }}}}_{\stackrel{\mathrm{‾}}{h}}$

Algorithm 4 Coarse-to-fine image resolution framework for inferring WV φ and attenuated backscatter χ.

Require:

(1) WV ${\widehat{\mathit{\phi }}}^{\text{init}}$ and atten. backs. ${\widehat{\stackrel{\mathrm{̃}}{\mathit{\chi }}}}^{\text{init}}$ initial values, (2) WV λ_w and atten. backs. λ_a tuning parameters and training photon counts ${\mathbf{Y}}_{\text{trn}}^{\left(\mathrm{on}\right)}$ and ${\mathbf{Y}}_{\text{trn}}^{\left(\mathrm{off}\right)}$, (3) coarsest resolution downsampling factor $\stackrel{\mathrm{‾}}{h}$.

1:

${\mathit{\phi }}^{(\stackrel{\mathrm{‾}}{h}+\mathrm{2})}={\widehat{\mathit{\phi }}}^{\text{init}},{\stackrel{\mathrm{̃}}{\mathit{\chi }}}^{(\stackrel{\mathrm{‾}}{h}+\mathrm{2})}={\widehat{\stackrel{\mathrm{̃}}{\mathit{\chi }}}}^{\text{init}}${Set initial values}

2:

$h=\stackrel{\mathrm{‾}}{h}$ {Iterate over the coarse-to-fine downsampling factors}

3:

while h≥1 do

4:

{Set initial values for current image resolution; downsample the initial WV estimate to the image resolution at which inference is done}

5:

$t=\mathrm{0},{\mathit{\varphi }}^{\left(t\right)}={\mathcal{D}}_h^{\downarrow }{\mathit{\phi }}^{(h+\mathrm{2})},{\mathit{\psi }}^{\left(t\right)}={\stackrel{\mathrm{̃}}{\mathit{\chi }}}^{(h+\mathrm{2})}$

6:

{Find minimizer of objective function (19) by alternating the minimization between the WV ϕ and atten. backs. ψ.}

7:

repeat

8:

{Solve for WV via adapted SPIRAL by minimizing Eq. (19); 𝒲 is the set of non-negative numbers}

$$\begin{array}{}\text{(28)}& \begin{array}{rl}{\mathit{\varphi }}^{(t+\mathrm{1})}=\arg \underset{\mathit{\varphi }\in \mathcal{W}}{min}{F}_{p_{\text{trn}},{\mathit{\lambda }}_{\mathrm{w}},{\mathit{\lambda }}_{\mathrm{a}}}(& {\mathcal{D}}_h^{\uparrow }\mathit{\varphi },{\mathit{\psi }}^{\left(t\right)};\\ & {\mathbf{Y}}_{\text{trn}}^{\left(\mathrm{on}\right)},{\mathbf{Y}}_{\text{trn}}^{\left(\mathrm{off}\right)})\end{array}\end{array}$$

9:

{Solve for atten. backs. via adapted SPIRAL by minimizing Eq. (19)}

$$\begin{array}{}\text{(29)}& \begin{array}{rl}{\mathit{\psi }}^{(t+\mathrm{1})}=\arg \underset{\mathit{\psi }}{min}{F}_{p_{\text{trn}},{\mathit{\lambda }}_{\mathrm{w}},{\mathit{\lambda }}_{\mathrm{a}}}(& {\mathcal{D}}_h^{\uparrow }{\mathit{\varphi }}^{(t+\mathrm{1})},\mathit{\psi };\\ & {\mathbf{Y}}_{\text{trn}}^{\left(\mathrm{on}\right)},{\mathbf{Y}}_{\text{trn}}^{\left(\mathrm{off}\right)})\end{array}\end{array}$$

10:

$t=t+\mathrm{1}$

11:

until stopping criteria is Eq. (24) met

12:

{Set initial values for finer image resolution}

13:

${\mathit{\phi }}^{\left(h\right)}={\mathcal{D}}_h^{\uparrow }{\mathit{\varphi }}^{\left(t\right)}$, ${\stackrel{\mathrm{̃}}{\mathit{\chi }}}^{\left(h\right)}={\mathit{\psi }}^{\left(t\right)}$, $h=max(h-\mathrm{2},\mathrm{1})$

14:

end while

15:

{Set final estimates of WV and atten. backs.}

16:

${\widehat{\mathit{\phi }}}_{\stackrel{\mathrm{‾}}{h},{\mathit{\lambda }}_{\mathrm{w}},{\mathit{\lambda }}_{\mathrm{a}}}={\mathit{\phi }}^{\left(\mathrm{1}\right)}$, ${\widehat{\stackrel{\mathrm{̃}}{\mathit{\chi }}}}_{\stackrel{\mathrm{‾}}{h},{\mathit{\lambda }}_{\mathrm{w}},{\mathit{\lambda }}_{\mathrm{a}}}={\stackrel{\mathrm{̃}}{\mathit{\chi }}}^{\left(\mathrm{1}\right)}$

17:

return  ${\widehat{\mathit{\phi }}}_{\stackrel{\mathrm{‾}}{h},{\mathit{\lambda }}_{\mathrm{w}},{\mathit{\lambda }}_{\mathrm{a}}}$, ${\widehat{\stackrel{\mathrm{̃}}{\mathit{\chi }}}}_{\stackrel{\mathrm{‾}}{h},{\mathit{\lambda }}_{\mathrm{w}},{\mathit{\lambda }}_{\mathrm{a}}}$

6 Experiment results

In this section we quantify the accuracy of the PTV-MPD WV measurements juxtaposed with the standard method; we use the mnemonics STND and PTV to make a distinction between the standard and PTV methods. We use RS WV measurements as an independent reference to estimate WV accuracies. The MPD with coincident RS observations was stationed (1) in proximity of the Atmospheric Radiation Measurement (ARM) SGP atmospheric observatory during spring and summer of 2019² and (2) at the NCAR Marshall field site in Boulder, Colorado, during winter and spring of 2020 and 2021 (Spuler et al., 2021). During the SGP MPD deployment, the ARM Raman lidar was making continuous WV measurements; Table 3 compares the laser power and telescope diameters of the Raman lidar and MPD.

Table 3 The first column indicates the lidar instrument and where the instrument was located. SGP refers to the Southern Great Plains in Oklahoma, and Marshall refers to Marshall field in Boulder, Colorado. The second and third column show the corresponding laser power and telescope diameter.

6.1 Radiosonde comparison methodology

Although the MPD instrument was in close proximity of the ARM site during the SGP deployment, the MPD, Raman lidar, and RS instruments all observe different atmospheric volumes at different altitudes, especially in highly convective conditions. Furthermore, the Raman lidar uses some RS observations for calibration, and therefore the instrument is not entirely independent of the RS (Newsom and Sivaraman, 2018). Therefore, we indicate the altitude at which the RS is horizontally more than 5 km away from the MPD, and we also show the ARM Raman lidar WV retrievals during the SGP deployment. In contrast to the Raman lidar, the MPD retrievals or calibrations do not employ RS observations. Instead, the MPD retrievals depend on approximate temperature and pressure profiles that are computed via an assumed lapse rate, and this rate is computed from a surface observation (Spuler et al., 2015; Nehrir et al., 2009).

6.2 Data selection methodology

For the experiments we created datasets that span over specific months given (1) the high variability of WV across different months and (2) that the MPD WV measurement capability is dependent on the low-altitude WV two-way transmittance. For all the datasets presented here, the observed photon counts are binned at range and time bins of 37.5 m and 5 min, respectively. The analyzed data from SGP are not comprehensive, and we instead selected a variety of challenging cases in an effort to identify potential issues in the PTV-MPD algorithm. Specifically, much of the data targeted instances where clouds and precipitation created challenging scenes to processes.

Specific to the MPD at SGP, the WV measurements start at 500 m range due to the poor geometric overlap below this range. The next-generation MPD deployed at the Marshall field site employed optically combined telescopes with narrow and wide fields of view to improve low-altitude geometric overlap, and as a result the WV measurements start at 150 m range (Spuler et al., 2021).

Table 4 The first column shows the names of the datasets used in the experiments. The second column lists the specific dates that are included in each dataset. The third and fourth columns list the number of days in each dataset and the total number of coincident RS profiles.

The first column of Table 4 lists the datasets, where the name of each dataset encodes the location with the year and month interval that each dataset covers. The second column lists the specific dates that are included in each dataset. The third and fourth columns list the number of days in each dataset and the total number of coincident RS profiles.

6.3 PTV-MPD tuning parameters and the coarse-to-fine (PTV-CF) configuration

The PTV-MPD method has the following input parameters, as indicated in Alg. 2:

1. the initial WV value ${\widehat{\mathit{\phi }}}^{\text{init}}$;

2. the WV λ_w and attenuated backscatter λ_a tuning parameters that set the degree to which the 2D piecewise constant regularization is promoted for the WV $\widehat{\mathit{\phi }}$ and attenuated backscatter $\widehat{\mathit{\chi }}$ geophysical variables (see Eq. 12);

3. the coarsest-resolution downsampling factor $\stackrel{\mathrm{‾}}{h}$, which controls at what image resolution the WV inference starts at.

In addition, the PTV-MPD method depends on an input cloud mask, which in itself is set by photon rate threshold (see Sect. 3.2). As discussed in Sect. 5.3.2, the initial WV value is set to zero [g m⁻³]. The WV λ_w and attenuated backscatter λ_a tuning parameters are determined through a cross-validation methodology (Friedman et al., 2001, chap. 7)(Oh et al., 2013); see Sect. 5.3.4 for more detail. For all of the experiments PTV-MPD uses the default set of TV regularizer tuning parameters, which are set between lines 1 to 5 in Alg. 2.

In regards to the resolution downsampling factor $\stackrel{\mathrm{‾}}{h}$, we juxtapose the PTV-MPD method with ($\stackrel{\mathrm{‾}}{h}>\mathrm{1}$) and without ($\stackrel{\mathrm{‾}}{h}=\mathrm{1}$) the coarse-to-fine configuration to test the hypothesis of whether the PTV-MPD coarse-to-fine image resolution framework does improve the WV due to inaccuracies in the initial attenuated backscatter value. The test error (22) is used to quantify the hypotheses testing, as indicated in Eq. (23). We use mnemonics, as indicated in Table 5, to make a distinction between the two PTV-MPD configurations and the standard method. The coarse-to-fine PTV-MPD method, PTV-CF, starts at a coarse image resolution, which is 9 times ($\stackrel{\mathrm{‾}}{h}=\mathrm{9}$) coarser than the finest image resolution.

Table 5 The mnemonics used to indicate the algorithm that was used to produce a water vapor estimate. The corresponding test errors are also shown.

In the following subsection we show individual PTV-MPD WV measurement results, and thereafter we present the WV measurement error statistics.

6.4 Individual retrieval results

Figures 3 and 4 show the WV measurements obtained using the standard and PTV-MPD methods (see Table 5) when the MPD was at SGP on 11 and 10 June 2019. Figure 5 show the WV measurements when the MPD was at Marshall on 10 December 2020. These examples have been selected to highlight instances of improvement and ongoing challenges for the PTV-CF and PTV methods. For all of these figures, panel (a) shows the attenuated backscatter χ. Panels (b) to (d) show the WV φ measurements of the STND, PTV-CF, and PTV methods. The rest of the panels, i.e., (e) and beyond, show the WV profiles of the different methods compared with that of the RS and Raman lidar; the vertical dashed white lines in panels (b) to (d) show the launch times of the RSs. Panel (a) is not masked to show what atmospheric features have been masked out in panels (b) to (d) by the gray areas.

Figure 3 The (a) attenuated-backscatter cross section and the water vapor (WV) measurements of the (b) standard method (STND) and the PTV-MPD method (c) with (PTV-CF) and (d) without (PTV) the coarse-to-fine configuration. Panel (a) is not masked to show what atmospheric features have been masked out in (b) to (d) by the gray areas. The dashed white vertical lines in (b) to (d) show the launch times of the radiosondes (RSs), where (e) to (h) show the RS and Raman lidar WV profiles. Panels (e) and (g) compare the course-to-fine PTV-MPD WV measurements against those of the RS and the standard method, whereas (f) and (h) show the same comparison for the PTV-MPD method without the course-to-fine configuration. The horizontal dashed lines in (e) to (h) show the altitude at which the RS was horizontally 5 km away from the MPD. The dashed boxes in (e) and (g) highlight regions that are discussed in the text.

Figure 4 The (a) attenuated backscatter cross section and the water vapor (WV) measurements of the (b) standard method (STND) and the PTV-MPD method (c) with (PTV-CF) and (d) without (PTV) the coarse-to-fine configuration. Panel (a) is not masked to show what atmospheric features have been masked out in (b) to (d) by the gray areas. The dashed vertical white lines in (b) to (d) show the launch times of the radiosondes (RSs), where (e) to (h) show the RS and Raman lidar WV profiles. Panels (e) and (g) compare the course-to-fine PTV-MPD WV measurements against those of the RS and the standard method, whereas (f) and (h) show the same comparison for the PTV-MPD method without the course-to-fine configuration. The horizontal dashed lines in (e) to (h) show the altitude at which the RS was horizontally 5 km away from the MPD. The dashed boxes in (e) and (g) highlight regions that are discussed in the text.

Figure 5 The (a) attenuated backscatter cross section and the water vapor (WV) measurements of the (b) standard method (STND) and the PTV-MPD method (c) with (PTV-CF) and (d) without (PTV) the coarse-to-fine configuration. Panel (a) is not masked to show what atmospheric features have been masked out in (b) to (d) by the gray areas. The dashed vertical white lines in (b) to (d) show the launch times of the radiosondes (RSs), where (e) and (f) show the RS WV profiles. Panel (e) compares the course-to-fine PTV-MPD WV measurements against those of the RS and the standard method, whereas (f) show the same comparison for the PTV-MPD method without the course-to-fine configuration. The horizontal dashed lines in (e) to (f) show the altitude at which the RS was horizontally 5 km away from the MPD. The dashed box in (e) highlights a region that is discussed in the text.

6.4.1 Comparing PTV and PTV-CF with STND

When comparing the PTV and PTV-CF WV measurements with that of the STND method, the STND method exhibits higher residual noise, particularly at higher altitudes, as indicated by Figs. 3i, ii, 4i, and 5i. These oscillatory residual noise artifacts in the STND WV measurements are due to the low-pass filter that is designed for high-fidelity WV measurements up to 4 km. In contrast, the PTV-CF and PTV methods produce more accurate WV measurements at higher altitudes, and these measurements do not exhibit large residual noise artifacts. The reasons why PTV-CF and PTV make more accurate high-altitude WV measurements are as follows:

1. PTV-CF and PTV employ the Poisson noise model with the MPD forward models to fit WV estimates onto the noisy observations, and

2. the PTV-CF and PTV methods, via the TV regularization, exploit the spatial and temporal correlations of the WV to separate the WV from the random noise.

In some cases, residual noise in the STND method may falsely appear to capture RS observed WV structure. An example of this is shown in Fig. 4ii, where an elevated WV layer appears to be better captured by the STND method than PTV-CF, where PTV-CF identified an elevated layer that consists of lower density and greater vertical extent than the RS and Raman lidar. The uncertainty analysis available with the standard method (i.e., STND), however, shows that the noise standard deviation in the WV estimate is greater than 2 g m⁻³, indicating that the apparent structure is not statistically significant and is predominantly noise. This fact also becomes more apparent when viewing the time-resolved plot in Fig. 4b.

Figure 4ii demonstrates the challenge in validating the MPD profiles, particularly for regions exhibiting significant structure at higher altitudes. In this instance, the Raman lidar starts to struggle measuring WV at 3.5 km altitude during the daytime, due to the high solar background radiation noise. In order to validate such cases, the MPD, Raman lidar, and RS instruments should ideally make WV measurements of the same atmospheric volume, which is not always possible. For the purposes of the statistical analysis presented in the next section, we will nevertheless treat the RS observations as “truth” despite possible discrepancies resulting from difference in atmospheric volume.

6.4.2 Comparing PTV-CF with PTV

Comparing Fig. 3c and d we see that the PTV WV measurement appears to be larger than that of PTV-CF between 09:00 and 12:00 UTC at 3 km. Figure 3iii quantitatively shows the large PTV WV measurement, which appears to be an artifact. In contrast, the corresponding PTV-CF WV measurement in Fig. 3e is closer to RS WV. The likely cause of the artifact produced by PTV is inaccuracies in the initial attenuated backscatter that are induced by

1. the assumed WV value when computing the initial value of the attenuated backscatter (see Sect. 5.3.2) and

2. long laser pulse, denoted by operator A in Eq. (5), that is not deconvolved from initial attenuated backscatter.

Figure 4ii and iii show another distinct WV measurement difference between the PTV-CF and PTV methods. The tenuous WV identified by PTV-CF correlates more with the RS WV measurement in magnitude and shape compared to the PTV method.

The higher accuracy of the PTV-CF WV measurements corresponding to Figs. 3iii and 4iii is because the initial values of the attenuated backscatter and WV are systematically improved with the coarse-to-fine image resolution framework. Specifically, while inferring the attenuated backscatter and WV from both the online and offline channels simultaneously at the coarsest image resolution, the PTV-CF method

1. implicitly deconvolves the long laser pulse from the initial attenuated backscatter

2. and constrains the inference of the WV to be at a coarse image resolution, which implicitly increases the SNR of the observations.

6.5 Water vapor measurement error statistics

The RS WV measurements are used as a reference to produce WV measurement error statistics of the different methods discussed in the following subsection. Next, we discuss the test errors between the PTV-CF and PTV methods. We then summarize the results of the error statistics.

6.5.1 Radiosonde comparisons

For each dataset we compute the root-mean-squared “error” (RMSE) per range and the relative RMSE per range; all the methods use the same mask to exclude the same WV pixels when computing the RMSE. The RMSE and relative RMSE were computed as follows. Let ${\widehat{\mathit{\phi }}}_n^{\left(\mathrm{1}\right)},{\widehat{\mathit{\phi }}}_n^{\left(\mathrm{2}\right)},\mathrm{\dots },{\widehat{\mathit{\phi }}}_n^{\left(L\right)}$ denote the WV profiles inferred by one of the methods, and let ${\mathit{\phi }}_n^{\left(\mathrm{1}\right)},{\mathit{\phi }}_n^{\left(\mathrm{2}\right)},\mathrm{\dots },{\mathit{\phi }}_n^{\left(L\right)}$ denote coincident RS WV profiles. The common mask profiles are denoted by ${\stackrel{\mathrm{̃}}{\mathbf{M}}}_n^{\left(\mathrm{1}\right)},{\stackrel{\mathrm{̃}}{\mathbf{M}}}_n^{\left(\mathrm{2}\right)},\mathrm{\dots },{\stackrel{\mathrm{̃}}{\mathbf{M}}}_n^{\left(L\right)}$; for example, ${\stackrel{\mathrm{̃}}{\mathbf{M}}}_n^{\left(\mathrm{1}\right)}=\mathrm{0}$ indicates that the WV measurement at row index n is masked out by either the standard or PTV-MPD method. The RMSE and relative RMSE (RRMSE) per dataset are computed by

$$\begin{array}{}\text{(30)}& {\displaystyle }{N}^{\left(l\right)}=\sum _{n=\mathrm{1}}^N{\stackrel{\mathrm{̃}}{\mathbf{M}}}_n^{\left(l\right)}\text{(31)}& {\displaystyle }\mathrm{RMSE}={\left({\displaystyle \frac{\mathrm{1}}{L}}\sum _{l=\mathrm{1},n=\mathrm{1}}^{L,N}{\displaystyle \frac{{\stackrel{\mathrm{̃}}{\mathbf{M}}}_n^{\left(l\right)}}{N^{\left(l\right)}}}{\left({\widehat{\mathit{\phi }}}_n^{\left(l\right)}-{\mathit{\phi }}_n^{\left(l\right)}\right)}^{\mathrm{2}}\right)}^{\mathrm{1}/\mathrm{2}},\end{array}$$

$$\begin{array}{}\text{(32)}& \begin{array}{rl}\mathrm{RRMSE}& =\mathrm{100}\times \mathrm{RMSE}\\ & \times {\left({\displaystyle \frac{\mathrm{1}}{L}}\sum _{l=\mathrm{1},n=\mathrm{1}}^{L,N}{\displaystyle \frac{{\stackrel{\mathrm{̃}}{\mathbf{M}}}_n^{\left(l\right)}}{N^{\left(l\right)}}}{\left({\mathit{\phi }}_n^{\left(l\right)}\right)}^{\mathrm{2}}\right)}^{-\mathrm{1}/\mathrm{2}}.\end{array}\end{array}$$

Figure 6 For the all the SGP datasets, (a) shows the root-mean-square “error” (RMSE) relative to the radiosonde (RS) water vapor (WV) measurements for the standard method (STND) and the PTV-MPD method with (PTV-CF) and without (PTV) the coarse-to-fine configuration, and (e) shows the corresponding relative RMSE. Panels (b) to (d) show the RMSE for each individual SGP dataset, and (f) to (h) show the corresponding relative RMSE.

Figure 7 For the all the Marshall datasets, (a) shows the root-mean-square “error” (RMSE) relative to the radiosonde (RS) water vapor (WV) measurements for the standard method (STND) and the PTV-MPD method with (PTV-CF) and without (PTV) the coarse-to-fine configuration, and (e) shows the corresponding relative RMSE. Panels (b) to (d) show the RMSE for each individual Marshall dataset, and (f) to (h) show the corresponding relative RMSE. Note that the scale of the horizontal axes in (a) to (d) are smaller than those shown in Fig. 6.

Figures 6a and 7a show the RMSE over the whole SGP and Marshall datasets, respectively, whereas Figs. 6e and 7e show the corresponding RRMSE. Figures 6b–d and 7b–d show the RMSE for each SGP and Marshall dataset (see Table 4), whereas Figs. 6f–h and 7f–h show the corresponding RRMSE.

From Figs. 6 and 7, we see that above 1.5 km the PTV-CF and PTV RMSEs and RRMSEs are comparable. Except for the lowest altitudes in Fig. 6c and d, the PTV-CF RMSEs are approximately less than 1 g m⁻³ for all altitudes, whereas the STND RMSE degrades with increasing altitude. For the SGP datasets, the PTV-CF and PTV RRMSEs are less than 100 % in most of the profiles up to 8 km altitude, where for the SGP 2019/04 dataset in Fig. 6f the RRMSE peaks at 125 % near 6 km altitude. The PTV-CF and PTV RRMSEs for the Marshall datasets in Fig. 7 are similar to the SGP RRMSEs below 8 km altitude.

The Raman lidar outperforms the PTV-CF and PTV when its observations have sufficiently high SNR. Otherwise, in high solar background conditions (i.e., lower SNR) the Raman lidar's effective range is diminished. By employing PTV-CF or PTV, the MPD, which had a power-aperture product approximately 500 times lower than the Raman lidar, is able to obtain higher-accuracy WV measurements above 4 km compared to the Raman lidar.

From Figs. 6a and 7a we see that the lowest-altitude WV measurements of the PTV-CF method are more accurate on average compared to the PTV and STND methods. These lower WV accuracies of the PTV and STND methods are due to lower near-range SNR observations and inaccuracies in the initial attenuated backscatter. As explained in Sect. 6.4.2, the course-to-fine methodology employed by PTV-CF mitigates the inaccuracies of the initial attenuated backscatter.

6.5.2 Test error comparisons

Table 6 shows that for five out of the six datasets the PTV-CF test errors $\left({\mathrm{err}}_{\mathrm{tst},\stackrel{\mathrm{‾}}{h}=\mathrm{9}}\right)$ are less than the PTV test error $\left({\mathrm{err}}_{\mathrm{tst},\stackrel{\mathrm{‾}}{h}=\mathrm{1}}\right)$. The lower test errors of PTV-CF correlate with the higher WV accuracy that PTV-CF is able to achieve compared to PTV at the lower altitudes.

Table 6 This table shows for which datasets the test error of the PTV-MPD method with the coarse-to-fine configuration $\left({\mathrm{err}}_{\mathrm{tst},\stackrel{\mathrm{‾}}{h}=\mathrm{9}}\right)$ is smaller (i.e., better) than PTV-MPD without the coarse-to-fine configuration $\left({\mathrm{err}}_{\mathrm{tst},\stackrel{\mathrm{‾}}{h}=\mathrm{1}}\right)$. The first two rows of the table, before the horizontal divider line, correspond to the total test error comparison over all the SGP and Marshall datasets. The rest of the rows correspond to the test error comparison over the specific SGP and Marshall datasets.

7 Conclusion and future work

PTV-MPD employs the photon-counting noise model to quantitatively fit the MPD forward models on the noisy observations, while encouraging accurate spatial and temporal correlations across all WV estimate pixels via the total variation regularizer function. This holistic approach allows for accurately measuring highly structured and varying WV fields at different altitudes and varying SNRs. In comparison, the standard method employs low-pass filters to reduce the residual noise in the WV estimate, and the low-pass filter bandwidth is optimized for lower-altitude (≤ 4 km) WV measurements. Therefore, the standard method, which uses high-altitude WV measurements, contains residual noise that degrades the measurements. By applying PTV-MPD to WV retrievals from the MPD, we have been able to extend the maximum altitude of the retrieval from 6 to 8 km, enabling a maximum coverage (after hardware improvements described in Spuler et al., 2021) from 150 m to 8 km by the instrument. In addition, the WV retrieval accuracy is substantially increased above 2 km. It is also notable that by employing the PTV-MPD method, the MPD WV measurement range extends beyond the that of the ARM-SGP Raman lidar, which has nearly 500 times the power-aperture product of the MPD (Newsom and Sivaraman, 2018).

We also demonstrated that without careful consideration of how PTV is adapted for the MPD instrument, low-altitude biases can be introduced in the PTV-MPD WV measurements. PTV-MPD requires an initial value of the attenuated backscatter, and any inaccuracies can induce biases in the PTV-MPD WV estimates. PTV-MPD can be made more robust against such biases by inferring the WV via a coarse-to-fine image resolution framework. By inferring the attenuated backscatter and WV at a coarse image resolution, inaccuracies in the initial attenuated backscatter are reduced, and subsequent finer image resolution attenuated backscatter estimates are more accurate, which allows for more accurate WV estimates.

As of now, PTV-MPD is computationally expensive since inferring the WV requires estimating the WV with several tuning parameters, and the optimal tuning parameter is selected through a cross-validation methodology. For example, with 144 CPU cores it takes 1 to 2 h to infer the WV using the PTV-MPD coarse-to-fine image resolution framework for 24 h of data; each CPU core estimates the WV for a specific tuning parameter. We are working towards developing a methodology to infer the optimal tuning parameter from the photon-counting observations, which will reduce the number of required CPU cores to 12 or less. In addition, adapting the PTV-MPD code to use a GPU instead of a CPU might reduce the computational time by at least 10-fold (Lee and Wright, 2008).

An additional benefit of decreasing the computational demands of PTV-MPD is that it would be able to process multiple days of consecutive observations and not just 24 h scenes as demonstrated in the results. A current workaround to process multiple-day scenes is to use a horizontal sliding window when inferring the WV. Specifically, the WV is inferred for consecutive overlapping 24 h periods, and the final multiple day WV image is obtained by averaging together the overlapping 24 h period WV estimates.

Future work includes the following goals.

1. The first is working towards quantifying the uncertainties of the PTV-MPD WV measurements using a bootstrapping methodology (Friedman et al., 2001). The uncertainty quantification will quantify how PTV-MPD WV measurements behave at the edges of scenes where less spatial information is available.

2. The second is investigating how PTV-MPD can be made more robust against saturated photon counts.

Appendix A: Jacobian matrix of MPD forward models (Eq. 5) and gradient matrix of loss function (Eq. 12)

Define the indicator function as

$$\begin{array}{}\text{(A1)}& \mathbb{I}\mathit{\{}k=l\mathit{\}}=\left\{\begin{array}{ll}\mathrm{1}& \text{ if }k=l\\ \mathrm{0}& \text{ if }k\ne l.\end{array}\right.\end{array}$$

The MPD forward model derivatives relative to WV and attenuated backscatter are

$$\begin{array}{}\text{(A2)}& {\displaystyle \frac{\partial }{\partial {\mathit{\phi }}_{j,l}}}{\mathbf{S}}_{n,k}^{\left(\mathit{\iota }\right)}(\mathit{\phi },\exp (\stackrel{\mathrm{̃}}{\mathit{\chi }}\left)\right)=-\mathrm{2}\mathrm{\Delta }r\mathrm{\Delta }t{\mathit{\sigma }}_{j,l}^{\left(\mathit{\iota }\right)}{\mathit{U}}_k\mathbb{I}\mathit{\{}k=l\mathit{\}},\text{(A3)}& {\displaystyle }\times \sum _{n^{\prime }=n}^{n+\mathrm{\Delta }N}{\mathbf{A}}_{n^{\prime },n}^{\mathrm{T}}\mathbb{I}\mathit{\{}j\le n^{\prime }\mathit{\}}{\stackrel{\mathrm{̃}}{\mathbf{S}}}_{n^{\prime },k}^{\left(\mathit{\iota }\right)}(\mathit{\phi },\exp (\stackrel{\mathrm{̃}}{\mathit{\chi }}\left)\right),\text{(A4)}& {\displaystyle \frac{\partial }{\partial {\stackrel{\mathrm{̃}}{\mathit{\chi }}}_{j,l}}}{\mathbf{S}}_{n,k}^{\left(\mathit{\iota }\right)}(\mathit{\phi },\exp (\stackrel{\mathrm{̃}}{\mathit{\chi }}\left)\right)=\mathrm{\Delta }t{\mathit{U}}_k\mathbb{I}\mathit{\{}k=l\mathit{\}},\text{(A5)}& {\displaystyle }\times \sum _{n^{\prime }=n}^{n+\mathrm{\Delta }N}{\mathbf{A}}_{n^{\prime },n}^{\mathrm{T}}\mathbb{I}\mathit{\{}j\le n^{\prime }\mathit{\}}{\stackrel{\mathrm{̃}}{\mathbf{S}}}_{n^{\prime },k}^{\left(\mathit{\iota }\right)}(\mathit{\phi },\exp (\stackrel{\mathrm{̃}}{\mathit{\chi }}\left)\right).\end{array}$$

The P-NLL loss function derivatives relative to WV and attenuated backscatter are

$$\begin{array}{}\text{(A6)}& {\displaystyle \frac{\partial }{\partial {\mathit{\phi }}_{j,l}}}{\mathcal{L}}_p\left({\mathbf{S}}^{\left(\mathit{\iota }\right)}\left(\mathit{\phi },\exp \left(\stackrel{\mathrm{̃}}{\mathit{\chi }}\right)\right);{\mathbf{Y}}^{\left(\mathit{\iota }\right)}\right)=-\mathrm{2}\mathrm{\Delta }r\mathrm{\Delta }tp{\mathit{\sigma }}_{j,l}^{\left(\mathit{\iota }\right)},\text{(A7)}& {\displaystyle }\times {\mathit{U}}_l\sum _{n^{\prime }=j}^{N+\mathrm{\Delta }N}{\stackrel{\mathrm{̃}}{\mathbf{S}}}_{n^{\prime },l}^{\left(\mathit{\iota }\right)}\left(\mathit{\phi },\exp \left(\stackrel{\mathrm{̃}}{\mathit{\chi }}\right)\right)\sum _{n=max(n^{\prime }-\mathrm{\Delta }N,\mathrm{1})}^{min(n^{\prime },N)}{\mathbf{A}}_{n^{\prime },n}^{\mathrm{T}},\text{(A8)}& {\displaystyle }\times {\mathit{M}}_{n,l}\left(\mathrm{1}-{\displaystyle \frac{{\mathbf{Y}}_{n,l}^{\left(\mathit{\iota }\right)}}{p{\mathbf{S}}_{n,l}^{\left(\mathit{\iota }\right)}\left(\mathit{\phi },\exp \left(\stackrel{\mathrm{̃}}{\mathit{\chi }}\right)\right)}}\right),\text{(A9)}& {\displaystyle \frac{\partial }{\partial {\stackrel{\mathrm{̃}}{\mathit{\chi }}}_{j,l}}}{\mathcal{L}}_p\left({\mathbf{S}}^{\left(\mathit{\iota }\right)}\left(\mathit{\phi },\exp \left(\stackrel{\mathrm{̃}}{\mathit{\chi }}\right)\right);{\mathbf{Y}}^{\left(\mathit{\iota }\right)}\right)=\mathrm{\Delta }tp,\text{(A10)}& {\displaystyle }\times {\mathit{U}}_l{\stackrel{\mathrm{̃}}{\mathbf{S}}}_{j,l}^{\left(\mathit{\iota }\right)}\left(\mathit{\phi },\exp \left(\stackrel{\mathrm{̃}}{\mathit{\chi }}\right)\right)\sum _{n=max(n^{\prime }-\mathrm{\Delta }N,\mathrm{1})}^{min(n^{\prime },N)}{\mathbf{A}}_{n^{\prime },n}^{\mathrm{T}},\text{(A11)}& {\displaystyle }\times {\mathit{M}}_{n,l}\left(\mathrm{1}-{\displaystyle \frac{{\mathbf{Y}}_{n,l}^{\left(\mathit{\iota }\right)}}{p{\mathbf{S}}_{n,l}^{\left(\mathit{\iota }\right)}\left(\mathit{\phi },\exp \left(\stackrel{\mathrm{̃}}{\mathit{\chi }}\right)\right)}}\right).\end{array}$$

Appendix B: Poisson negative log likelihood Eq. 8 properties

Here we show that for high solar background radiation (i.e., lower SNR) of the photon-counting images there can be multiple estimates of the WV φ or attenuated backscatter $\stackrel{\mathrm{̃}}{\mathit{\chi }}$ that minimize the objective function (14). Specifically, we show that the P-NLL (8) has more the one minimizer. Without loss of generality we assume here that

1. the laser pulse duration is equal to the sampling intervals, meaning that ΔN=0 and the matrix A is just an identity matrix,

2. none of the photon counts are masked out by matrix M,

3. and p=1.

With these assumptions the Hessian matrix for observation column k of the P-NLL (8) relative to WV φ for given attenuated backscatter χ equals

$$\begin{array}{}\text{(B1)}& {\displaystyle }{\mathit{H}}_{j,l}(\mathit{\phi };\mathit{\chi })=\mathrm{4}\mathrm{\Delta }{r}^{\mathrm{2}}{\mathit{\sigma }}_{j,k}^{\left(\mathit{\iota }\right)}{\mathit{\sigma }}_{l,k}^{\left(\mathit{\iota }\right)}\sum _{l^{\prime }=l}^N\mathbb{I}\left\{j\le l^{\prime }\right\},\text{(B2)}& {\displaystyle }\times \left({\stackrel{\mathrm{̃}}{\mathit{S}}}_{l^{\prime },k}^{\left(\mathit{\iota }\right)}(\mathit{\phi },\mathit{\chi })-{\mathit{U}}_k\mathrm{\Delta }t{\mathit{b}}_k\right),\text{(B3)}& {\displaystyle }\times \left(\mathrm{1}-{\mathbf{Y}}_{l^{\prime },k}^{\left(\mathit{\iota }\right)}{\displaystyle \frac{{\mathit{U}}_k\mathrm{\Delta }t{\mathit{b}}_k}{{\left[{\stackrel{\mathrm{̃}}{\mathit{S}}}_{l^{\prime },k}^{\left(\mathit{\iota }\right)}(\mathit{\phi },\mathit{\chi })\right]}^{\mathrm{2}}}}\right).\end{array}$$

If this Hessian matrix is positive definite, meaning all its eigenvalues are positive, then it implies that for a fixed attenuated backscatter χ there is a unique WV φ estimate that minimizes the P-NLL and therefore the objective function (14) (Boyd and Vandenberghe, 2004). The Hessian matrix is not positive definite if

$$\begin{array}{}\text{(B4)}& {\stackrel{\mathrm{̃}}{\mathit{S}}}_{l^{\prime },k}^{\left(\mathit{\iota }\right)}(\mathit{\phi },\mathit{\chi })\le {\left[{\mathbf{Y}}_{l^{\prime },k}^{\left(\mathit{\iota }\right)}{\mathit{U}}_k\mathrm{\Delta }t{\mathit{b}}_k\right]}^{\frac{\mathrm{1}}{\mathrm{2}}}.\end{array}$$

Whether this inequality is met is dependent on the attenuated backscatter, WV transmittance, and laser energy. For example, for humid summer days the WV transmittance will significantly decrease the SSLE relative to the solar background radiation, consequently limiting PTV-MPD's ability to uniquely infer the WV at higher altitudes.

In the cases where the backscattered photons are comparable to the observed background counts according to Eq. (B4), we leave it to future work to determine the altitudes at which PTV-MPD can uniquely infer the WV.