the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Extending water vapor measurement capability of photonlimited differential absorption lidars through simultaneous denoising and inversion
Willem J. Marais
Matthew Hayman
The micropulse differential absorption lidar (MPD) was developed at Montana State University (MSU) and the National Center for Atmospheric Research (NCAR) to perform rangeresolved water vapor (WV) measurements using lowpower lasers and photoncounting detectors. The MPD has proven to produce accurate WV measurements up to 6 km altitude. However, the MPD's ability to produce accurate higheraltitude 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 groundbased 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 PTVMPD, 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. PTVMPD achieves the improvement by simultaneously denoising the MPD noisy observations and inferring the WV by separating the random noise from the nonrandom WV.
An analysis with 130 radiosonde (RS) comparisons shows that the relative rootmeansquare difference (RRMSE) of WV measurements between RS and PTVMPD 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 PTVMPD, 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 poweraperture product 500 times greater than that of the MPD.
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 rangeresolved 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 rangeresolved 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 lowpower, lowcost, highreliability diode lasers that enable class 1M eyesafe 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 firstprinciplesbased 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 higheraltitude 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 lowpass filters on the algebraically computed WV using a Gaussian smoothing kernel or Savitzky–Golay filter (Schafer, 2011). The signaltonoise ratio (SNR) of the WV measurements can be further improved upon by reducing the vertical and horizontal resolutions of the photoncounting observations. However, reducing the resolutions introduces systematic biases in the WV due to the nonlinearity of the singlescatter 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 oversmoothed (resulting in biases), while other parts are undersmoothed (resulting in random error) (Hayman et al., 2020).
Figure 1 illustrates the shortcomings of lowpass filtering, where Fig. 1b and e show the standard method estimated WV image and profile obtained using lowpass filtering. The profile is juxtaposed with a coincident radiosonde (RS) WV profile. Although the bandwidth of the lowpass 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 lowpass filter also oversmooths rapidly changing WV features such as embedded dry regions.
1.2 Proposed approach to extending MPD capability
Advances have been made in denoising photoncounting 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 photoncounting 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 groundbased photoncounting 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 PTVMPD, of the same scene and profile of Fig. 1c and e; using the range intervals indicated by Fig. 1i and ii, the PTVMPD WV measurements are more accurate at higher altitudes and at the embedded dry region compared to the standard method.
The PTVMPD 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 tradeoff of blurring the backscatter optical signal described by the standard singlescatter lidar equation (SSLE).

Simultaneous inference. PTV for the HSRL infers the backscatter and extinction cross sections with a twostep approach in order to isolate geometric overlap calibration biases to the extinction cross section. However, with PTVMPD a twostep 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 PTVMPD 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 twostep 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 PTVMPD is sensitive to nearrange inaccuracies in the initial attenuatedbackscatter estimate, which induces inaccuracies in the nearrange WV measurements. Therefore, the third adaptation involves making PTVMPD more robust against inaccuracies in the initial attenuated backscatter.
1.3 Contributions
The immediate contributions of this paper are the

adaptation of the PTV method (i.e., PTVMPD) for the DIAL technique and

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 PTVMPD methods are discussed in Sects. 4 and 5. The paper ends with experiment results in Sect. 6 and the conclusion thereafter.
Geophysical variables and the lidar forward models are written as matrices. When we introduce a geophysicsrelated 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.
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 midduration 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{\u0303}}{\mathit{r}}}_{n}$, where $n=\mathrm{1},\mathrm{2},\mathrm{\dots},N$ indexes the observation range axis; the relation between the observation ${\stackrel{\mathrm{\u0303}}{\mathit{r}}}_{n}$ and measurement ${\mathit{r}}_{{n}^{\prime}}$ range axes is
The SSLE is a function of the unknown WV ${\mathit{\phi}}_{{n}^{\prime},k}$ [g m^{−3}] and uncalibrated attenuated backscatter ${\mathit{\chi}}_{{n}^{\prime},k}$ [m^{−1} sr^{−1}], 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 singlescatter lidar matrix function at measurement range ${\mathit{r}}_{{n}^{\prime}}$ for channel ι is denoted by
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^{3}], and geometric overlap ${\mathit{O}}_{{n}^{\prime}}$. The WV absorption [m^{2} g^{−1}] of channel ι is denoted by ${\mathit{\sigma}}_{{n}^{\prime},k}^{\left(\mathit{\iota}\right)}$ and is precomputed 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
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{\u0303}}{\mathit{r}}}_{n}$ for channel ι is defined by
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 coaxial 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_{1}, 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_{1}. Hence, the WV φ_{1,k} at range r_{1} represents the rangeaverage WV from range 0 up to r_{1}.
3.1 Photoncounting noise probability mass function
The photoncounting observations at range ${\stackrel{\mathrm{\u0303}}{\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
where
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 Photoncounting noise negative log likelihood
The Poisson negative log likelihood (PNLL) is used by the PTVMPD method to quantify the fitting of the WV and attenuated backscatter onto the noisy observations via the MPD forward model Eq. (5). The PNLL function of channel ι is denoted by
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.
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.
The photon count profiles first undergo lowpass 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 lowpassfiltered 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 lowpass filter is fixed across the scene. Thus, the WV estimates are often oversmoothed in highly dynamic, highSNR regions and undersmoothed 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 attenuatedbackscatter 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.
The PTV method, originally developed for photoncounting 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

the photoncounting noise can be accurately modeled by the Poisson PMF in Eq. (6),

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

the geophysical variable (i.e., WV) image that we want to estimate can be accurately approximated with a twodimensional (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 photoncounting 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 PNLL 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) PNLL composited with the forward model and (2) the TV penalty functions. The PNLL 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
where X ≡ geophysical variable, Y ≡ photoncounting observations, $\stackrel{\mathrm{\u0303}}{\mathrm{\ell}}(\mathbf{X};\mathbf{Y})$ ≡ PNLL composited with forward model, e.g., Eq. (8).
The (anisotropic) TV penalty function is defined as (Harmany et al., 2012)
and is weighted with a tuning parameter $\stackrel{\mathrm{\u0303}}{\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 crossvalidation methodology (Friedman et al., 2001, chap. 7)(Oh et al., 2013), and the parameter is not set by “expert opinion”.
5.2 PTVMPD: 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 PTVMPD 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 attenuatedbackscatter 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 twostep 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 twostep 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 attenuatedbackscatter value. Specifically, when computing the initial value of the attenuated backscatter algebraically from the photoncounting 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 PTVMPD employ a coarsetofine image resolution framework when inferring the WV (Marais and Willett, 2017). With the coarsetofine 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 closerange observations are implicitly increased. In a sequence from coarsetofine image resolution, we use the coarseresolution WV and attenuatedbackscatter estimates as initial values when inferring the finer image resolution WV and attenuated backscatter.
5.3 PTVMPD method algorithmic details
Using the conceptual optimization framework in Eq. (10) as a reference, the PTVMPD 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 PTVMPD 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 coarsetofine 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 coarsetofine image resolution framework that addresses adaption (A3).
The formulation of the PTVMPD includes tuning parameters that are computed using a crossvalidation (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 coarsetofine 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 coarsetofine image resolution. The test error computation is also discussed in Sect. 5.3.4.
Figure 2 gives a broad pictorial overview of how each part of the PTVMPD framework is implemented with modularized subalgorithms; the purposes of the subalgorithms 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.
These algorithms are outlined in Sect. 5.3.5, which includes specific details about the SPIRAL method adaptations.
5.3.1 Preliminary for coarsetofine framework: the loss and objective functions
The PTVMPD loss function differs from the PTV loss function in two respects; cf. Eq. (23) in Marais et al. (2016). First, the loss function of PTVMPD is defined such that we infer the log _{e} of the attenuated backscatter denoted by ${\stackrel{\mathrm{\u0303}}{\mathit{\chi}}}_{n,k}\equiv {\mathrm{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 PTVMPD loss function as the sum of the normalized PNLL function of each channel, where the normalization ensures that each PNLL 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 PTVMPD loss function is defined by
which is the weighted sum of the PNLL of each channel where each PNLL 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
The objective function that is minimized by PTVMPD is defined as
which is parameterized with the WV φ and attenuated backscatter $\stackrel{\mathrm{\u0303}}{\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 nonnegative, 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 lowSNR regions of the photoncounting images there can be multiple estimates of the WV φ or attenuated backscatter $\stackrel{\mathrm{\u0303}}{\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 PTVMPD will require accurate initial values of the WV φ and attenuated backscatter $\stackrel{\mathrm{\u0303}}{\mathit{\chi}}$ whenever the photon count SNRs are dominated by the background rate b^{(ι)} such as at close range and in highaltitude 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^{−3}]; this WV initial value reflects that we assume no a priori information about the WV other than (1) the assumed lapse rate^{1} 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{\u0303}}{\mathit{\chi}}}}^{\text{init}}$, is computed from the offline channel with the initial WV ${\widehat{\mathit{\phi}}}^{\text{init}}$ via
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 coarsetofine inference framework
With the coarsetofine image resolution inference framework, we start with estimating a coarseresolution WV image and use the coarseresolution estimate as an initial WV value at a finer image resolution; this technique has proven useful to more accurately denoising and inverting images for lowSNR photonlimited 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 coarsetofine 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.

A series of values $\stackrel{\mathrm{\u203e}}{h}={h}_{\mathrm{1}}>{h}_{\mathrm{2}}>\mathrm{\dots}>{h}_{L}=\mathrm{1}$ that represents the downsampling factors of the coarsetofine image resolutions at which the WV will be inferred; $\stackrel{\mathrm{\u203e}}{h}$ represents the coarsest image resolution.

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 coarseresolution 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 coarsetofine image resolution index. The coarsetofine 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 attenuatedbackscatter estimate ${\mathit{\psi}}_{{h}_{l}}$. In addition, define the current image resolution objective function.
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):
where 𝒲 is the set of nonnegative numbers.
Final step. If l=L, the final WV and attenuatedbackscatter estimates are denoted by ${\widehat{\mathit{\phi}}}_{\stackrel{\mathrm{\u203e}}{h},{\mathit{\lambda}}_{\mathrm{w}},{\mathit{\lambda}}_{\mathrm{a}}}$ and ${\widehat{\stackrel{\mathrm{\u0303}}{\mathit{\chi}}}}_{\stackrel{\mathrm{\u203e}}{h},{\mathit{\lambda}}_{\mathrm{w}},{\mathit{\lambda}}_{\mathrm{a}}}$ meaning
and
The subscripts $\stackrel{\mathrm{\u203e}}{h}$, λ_{w} and λ_{a} indicate that the estimates of the WV φ and attenuated backscatter $\stackrel{\mathrm{\u0303}}{\mathit{\chi}}$ are the specific to coarsetofine image resolution configuration parameterized by $\stackrel{\mathrm{\u203e}}{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{\u0303}}{\mathit{\chi}}$ estimates against noisy observations that are statistically independent of these estimates. This is achieved by thinning the nonmasked 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 nonmasked pixel in a photoncounting image Y^{(ι)}, the Poisson thinning technique randomly samples the individual photon counts, and the sampled photon counts are placed in three photoncounting 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)}$ photoncounting 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{\u203e}}{h},{\mathit{\lambda}}_{\mathrm{w}},{\mathit{\lambda}}_{\mathrm{a}}}$ and attenuated backscatter ${\widehat{\stackrel{\mathrm{\u0303}}{\mathit{\chi}}}}_{\stackrel{\mathrm{\u203e}}{h},{\mathit{\lambda}}_{\mathrm{w}},{\mathit{\lambda}}_{\mathrm{a}}}$ for specific WV λ_{w} and attenuated backscatter λ_{a} tuning parameters and coarsetofine configuration $\stackrel{\mathrm{\u203e}}{h}$. Specifically, when using the training photon counts the objective function in Eq. (16) will be equal to
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 PNLL (8) and the forward model (5) and is denoted by
The optimum tuning parameters corresponds to the smallest validation error, and this give us the WV and attenuated backscatter estimates
Test photon counts
The test error for the estimates ${\widehat{\mathit{\phi}}}_{\stackrel{\mathrm{\u203e}}{h}}$ and ${\widehat{\stackrel{\mathrm{\u0303}}{\mathit{\chi}}}}_{\stackrel{\mathrm{\u203e}}{h}}$ are also computed with the PNLL (8) and the forward model (5) and is denoted by
To test the hypotheses as to whether the coarsetofine image resolution framework does improve the WV measurements, we evaluate the inequality
5.3.5 The PTVMPD algorithm
Algorithm 2 and its subalgorithms (Algs. 3 and 4) outline the PTVMPD algorithms as described by Fig. 2. Algorithm 2 takes as input the photon counts and the coarse image resolution at which PTVMPD will start inferring the WV; the starting image resolution is expressed as a downsampling factor $\stackrel{\mathrm{\u203e}}{h}\ge \mathrm{1}$. Algorithm 3 outlines the CV methodology, and Alg. 4 outlines the coarsetofine 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.,
is less than 10^{−5} then the loop in Alg. 4 terminates, where
In this section we quantify the accuracy of the PTVMPD 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^{2} 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.
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 lowaltitude WV twoway 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 PTVMPD 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 nextgeneration MPD deployed at the Marshall field site employed optically combined telescopes with narrow and wide fields of view to improve lowaltitude geometric overlap, and as a result the WV measurements start at 150 m range (Spuler et al., 2021).
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 PTVMPD tuning parameters and the coarsetofine (PTVCF
) configuration
The PTVMPD method has the following input parameters, as indicated in Alg. 2:

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

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);

the coarsestresolution downsampling factor $\stackrel{\mathrm{\u203e}}{h}$, which controls at what image resolution the WV inference starts at.
In addition, the PTVMPD 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^{−3}]. The WV λ_{w} and attenuated backscatter λ_{a} tuning parameters are determined through a crossvalidation methodology (Friedman et al., 2001, chap. 7)(Oh et al., 2013); see Sect. 5.3.4 for more detail. For all of the experiments PTVMPD 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{\u203e}}{h}$, we juxtapose the PTVMPD method with ($\stackrel{\mathrm{\u203e}}{h}>\mathrm{1}$) and without ($\stackrel{\mathrm{\u203e}}{h}=\mathrm{1}$) the coarsetofine configuration to test the hypothesis of whether the PTVMPD coarsetofine 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 PTVMPD configurations and the standard method. The coarsetofine PTVMPD method, PTVCF
, starts at a coarse image resolution, which is 9 times ($\stackrel{\mathrm{\u203e}}{h}=\mathrm{9}$) coarser than the finest image resolution.
In the following subsection we show individual PTVMPD 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 PTVMPD 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 PTVCF
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
, PTVCF
, 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.
6.4.1 Comparing PTV
and PTVCF
with STND
When comparing the PTV
and PTVCF
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 lowpass filter that is designed for highfidelity WV measurements up to 4 km. In contrast, the PTVCF
and PTV
methods produce more accurate WV measurements at higher altitudes, and these measurements do not exhibit large residual noise artifacts. The reasons why PTVCF
and PTV
make more accurate highaltitude WV measurements are as follows:

PTVCF
andPTV
employ the Poisson noise model with the MPD forward models to fit WV estimates onto the noisy observations, and 
the
PTVCF
andPTV
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 PTVCF
, where PTVCF
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^{−3}, indicating that the apparent structure is not statistically significant and is predominantly noise. This fact also becomes more apparent when viewing the timeresolved 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 PTVCF
with PTV
Comparing Fig. 3c and d we see that the PTV
WV measurement appears to be larger than that of PTVCF
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 PTVCF
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

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

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 PTVCF
and PTV
methods. The tenuous WV identified by PTVCF
correlates more with the RS WV measurement in magnitude and shape compared to the PTV
method.
The higher accuracy of the PTVCF
WV measurements corresponding to Figs. 3iii and 4iii is because the initial values of the attenuated backscatter and WV are systematically improved with the coarsetofine 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 PTVCF
method

implicitly deconvolves the long laser pulse from the initial attenuated backscatter

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 PTVCF
and PTV
methods. We then summarize the results of the error statistics.
6.5.1 Radiosonde comparisons
For each dataset we compute the rootmeansquared “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{\u0303}}{\mathbf{M}}}_{n}^{\left(\mathrm{1}\right)},{\stackrel{\mathrm{\u0303}}{\mathbf{M}}}_{n}^{\left(\mathrm{2}\right)},\mathrm{\dots},{\stackrel{\mathrm{\u0303}}{\mathbf{M}}}_{n}^{\left(L\right)}$; for example, ${\stackrel{\mathrm{\u0303}}{\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 PTVMPD method. The RMSE and relative RMSE (RRMSE) per dataset are computed by
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 PTVCF
and PTV
RMSEs and RRMSEs are comparable. Except for the lowest altitudes in Fig. 6c and d, the PTVCF
RMSEs are approximately less than 1 g m^{−3} for all altitudes, whereas the STND
RMSE degrades with increasing altitude. For the SGP datasets, the PTVCF
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 PTVCF
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 PTVCF
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 PTVCF
or PTV
, the MPD, which had a poweraperture product approximately 500 times lower than the Raman lidar, is able to obtain higheraccuracy WV measurements above 4 km compared to the Raman lidar.
From Figs. 6a and 7a we see that the lowestaltitude WV measurements of the PTVCF
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 nearrange SNR observations and inaccuracies in the initial attenuated backscatter. As explained in Sect. 6.4.2, the coursetofine methodology employed by PTVCF
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 PTVCF
test errors $\left({\mathrm{err}}_{\mathrm{tst},\stackrel{\mathrm{\u203e}}{h}=\mathrm{9}}\right)$ are less than the PTV
test error $\left({\mathrm{err}}_{\mathrm{tst},\stackrel{\mathrm{\u203e}}{h}=\mathrm{1}}\right)$. The lower test errors of PTVCF
correlate with the higher WV accuracy that PTVCF
is able to achieve compared to PTV
at the lower altitudes.
PTVMPD employs the photoncounting 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 lowpass filters to reduce the residual noise in the WV estimate, and the lowpass filter bandwidth is optimized for loweraltitude (≤ 4 km) WV measurements. Therefore, the standard method, which uses highaltitude WV measurements, contains residual noise that degrades the measurements. By applying PTVMPD 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 PTVMPD method, the MPD WV measurement range extends beyond the that of the ARMSGP Raman lidar, which has nearly 500 times the poweraperture product of the MPD (Newsom and Sivaraman, 2018).
We also demonstrated that without careful consideration of how PTV is adapted for the MPD instrument, lowaltitude biases can be introduced in the PTVMPD WV measurements. PTVMPD requires an initial value of the attenuated backscatter, and any inaccuracies can induce biases in the PTVMPD WV estimates. PTVMPD can be made more robust against such biases by inferring the WV via a coarsetofine 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, PTVMPD is computationally expensive since inferring the WV requires estimating the WV with several tuning parameters, and the optimal tuning parameter is selected through a crossvalidation methodology. For example, with 144 CPU cores it takes 1 to 2 h to infer the WV using the PTVMPD coarsetofine 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 photoncounting observations, which will reduce the number of required CPU cores to 12 or less. In addition, adapting the PTVMPD code to use a GPU instead of a CPU might reduce the computational time by at least 10fold (Lee and Wright, 2008).
An additional benefit of decreasing the computational demands of PTVMPD 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 multipleday 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.

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

The second is investigating how PTVMPD can be made more robust against saturated photon counts.
Define the indicator function as
The MPD forward model derivatives relative to WV and attenuated backscatter are
The PNLL loss function derivatives relative to WV and attenuated backscatter are
Here we show that for high solar background radiation (i.e., lower SNR) of the photoncounting images there can be multiple estimates of the WV φ or attenuated backscatter $\stackrel{\mathrm{\u0303}}{\mathit{\chi}}$ that minimize the objective function (14). Specifically, we show that the PNLL (8) has more the one minimizer. Without loss of generality we assume here that

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

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

and p=1.
With these assumptions the Hessian matrix for observation column k of the PNLL (8) relative to WV φ for given attenuated backscatter χ equals
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 PNLL and therefore the objective function (14) (Boyd and Vandenberghe, 2004). The Hessian matrix is not positive definite if
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 PTVMPD'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 PTVMPD can uniquely infer the WV.
The PTVMPD code will be made available through a Git repository https://github.com/WillemMarais/poissontotalvariation (last access: 26 August 2022) in the near future.
All of the MPD data products used in this work were produced by NCAR and are available upon request; the raw MPD data are available at https://doi.org/10.26023/MX0DZ722M406 (NCAR/EOL MPD Team, 2020). The SGP radiosonde data are available via https://doi.org/10.5439/1373934 (ARM, 2022). The Marshall radiosonde data are available via the GCOS Reference UpperAir Network (GRUAN) website https://www.gruan.org (last access: 26 August 2022, Hurst et al., 2022).
This work was a collaboration between the authors. WJM is the original developer of PTV and provided expertise in advanced signal processing. MH provided expertise in optical system modeling and DIAL observation. The PTVMPD development was the result of a collaborative effort to mutually leverage the two disciplines.
The contact author has declared that none of the authors has any competing interests.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
The authors would like to thank Robert A. Stillwell and Scott M. Spuler for their feedback on this paper. This material is based upon work supported by NSF AGS1930907 and the National Center for Atmospheric Research, which is a major facility sponsored by the National Science Foundation under cooperative agreement no. 1852977. We thank Dale Hurst, Patrick Cullis, Emrys Hall, and Allen Jordan of the NOAA Global Monitoring Laboratory for providing the sounding data from the Marshall field site.
This research has been supported by the Division of Atmospheric and Geospace Sciences (grant no. NSF AGS1930907) and the National Science Foundation (cooperative agreement no. 1852977). The MPD Net demo field campaign was supported by the Office of Biological and Environmental Research of the U.S. Department of Energy (under grant no. DOE/SCARM20002) as part of the Atmospheric Radiation Measurement (ARM) user facility, an Office of Science user facility.
This paper was edited by Hartwig Harder and reviewed by three anonymous referees.
Ahn, S. and Fessler, J.: Globally convergent image reconstruction for emission tomography using relaxed ordered subsets algorithms, IEEE T. Med. Imaging, 22, 613–626, https://doi.org/10.1109/TMI.2003.812251, 2003. a
ARM: Raman Lidar: Aerosol backscatter, scattering ratio, lidar ratio, extinction, cloud mask, and linear depolarization ratio derived from Thorson FEX code, ARM (Atmospheric Radiation Measurement) user facility [data set], https://doi.org/10.5439/1373934, 2022. a
Azzari, L. and Foi, A.: Variance stabilization in Poisson image deblurring, in: 2017 IEEE 14th International Symposium on Biomedical Imaging (ISBI 2017), 18–21 April 2017, Melbourne, VIC, Australia, 728–731, https://doi.org/10.1109/ISBI.2017.7950622, 2017. a, b
Beck, A. and Tetruashvili, L.: On the Convergence of Block Coordinate Descent Type Methods, SIAM J. Optimiz., 23, 2037–2060, https://doi.org/10.1137/120887679, 2013. a, b
Boyd, S. and Vandenberghe, L.: Convex optimization, Cambridge university press, ISBN 9780521833783, 2004. a
Cinlar, E.: Introduction to stochastic processes, Dover Publications, ISBN 9780486497976, 2013. a
Donovan, D. P., Whiteway, J. A., and Carswell, A. I.: Correction for nonlinear photoncounting effects in lidar systems, Appl. Opt., 32, 6742–6753, https://doi.org/10.1364/AO.32.006742, 1993. a, b
Fessler, J. A.: Optimization Methods for Magnetic Resonance Image Reconstruction: Key Models and Optimization Algorithms, IEEE Signal Proc. Mag., 37, 33–40, https://doi.org/10.1109/MSP.2019.2943645, 2020. a
Friedman, J., Hastie, T., and Tibshirani, R.: The elements of statistical learning, vol. 1, Springer series in statistics, New York, ISBN 0387848576, 2001. a, b, c
Geerts, B., Parsons, D., Ziegler, C. L., Weckwerth, T. M., Turner, D. D., Wurman, J., Kosiba, K., Rauber, R. M., McFarquhar, G. M., Parker, M. D., Schumacher, R. S., Coniglio, M. C., Haghi, K., Biggerstaff, M. I., Klein, P. M., Jr., W. A. G., Demoz, B. B., Knupp, K. R., Ferrare, R. A., Nehrir, A. R., Clark, R. D., Wang, X., Hanesiak, J. M., Pinto, J. O., and Moore, J. A.: The 2015 Plains Elevated Convection At Night (PECAN) field project, B. Am. Meteorol. Soc., 98, 767–786, https://doi.org/10.1175/BAMSD1500257.1, 2016. a
Harmany, Z. T., Marcia, R. F., and Willett, R. M.: This is SPIRALTAP: Sparse Poisson Intensity Reconstruction ALgorithms – Theory and Practice, IEEE T. Image Process., 21, 1084–1096, https://doi.org/10.1109/TIP.2011.2168410, 2012. a, b, c, d, e, f, g, h
Hayman, M. and Spuler, S.: Demonstration of a diodelaserbased high spectral resolution lidar (HSRL) for quantitative profiling of clouds and aerosols, Opt. Express, 25, A1096–A1110, https://doi.org/10.1364/OE.25.0A1096, 2017. a, b
Hayman, M., Stillwell, R. A., and Spuler, S. M.: Optimization of linear signal processing in photon counting lidar using Poisson thinning, Opt. Lett., 45, 5213, https://doi.org/10.1364/OL.396498, 2020. a, b, c
Hurst, D., Cullis, P., Hall, E., and Jordan, A.: GCOS Reference UpperAir Network (GRUAN), GRUAN Lead Centre, DWD, https://www.gruan.org, last access: 26 August 2022. a
Jensen, M. P., Petersen, W. A., Bansemer, A., Bharadwaj, N., Carey, L. D., Cecil, D. J., Collis, S. M., Genio, A. D. D., Dolan, B., Gerlach, J., Giangrande, S. E., Heymsfield, A., Heymsfield, G., Kollias, P., Lang, T. J., Nesbitt, S. W., Neumann, A., Poellot, M., Rutledge, S. A., Schwaller, M., Tokay, A., Williams, C. R., Wolff, D. B., Xie, S., and Zipser, E. J.: The Midlatitude Continental Convective Clouds Experiment (MC3E), B. Am. Meteorol. Soc., 97, 1667–1686, https://doi.org/10.1175/BAMSD1400228.1, 2016. a
Kelley, C. T.: Iterative Methods for Optimization, Society for Industrial and Applied Mathematics, https://doi.org/10.1137/1.9781611970920, 1999. a
Lee, S. and Wright, S. J.: Implementing algorithms for signal and image reconstruction on graphical processing units, Optimization Online, 10, University of WisconsinMadison, https://optimizationonline.org/2008/10/2131/ (last access: 2 September 2022), 2008. a
Marais, W. and Willett, R.: ProximalGradient methods for poisson image reconstruction with BM3DBased regularization, in: 2017 IEEE 7th International Workshop on Computational Advances in MultiSensor Adaptive Processing (CAMSAP), 10–13 December 2017, Curacao, 1–5, https://doi.org/10.1109/CAMSAP.2017.8313128, 2017. a, b, c
Marais, W. J., Holz, R. E., Hu, Y. H., Kuehn, R. E., Eloranta, E. E., and Willett, R. M.: Approach to simultaneously denoise and invert backscatter and extinction from photonlimited atmospheric lidar observations, Appl. Opt., 55, 8316–8334, https://doi.org/10.1364/AO.55.008316, 2016. a, b, c, d, e, f, g
Müller, J. W.: Deadtime problems, Nucl. Instrum. Meth., 112, 47–57, https://doi.org/10.1016/0029554X(73)907738, 1973. a
National Academies of Sciences, Engineering, and Medicine: The Future of Atmospheric Boundary Layer Observing, Understanding, and Modeling: Proceedings of a Workshop, The National Academies Press, Washington, DC, https://doi.org/10.17226/25138, 2018a. a
National Academies of Sciences, Engineering, and Medicine: Thriving on Our Changing Planet: A Decadal Strategy for Earth Observation from Space, The National Academies Press, Washington, DC, https://doi.org/10.17226/24938, 2018b. a
NCAR/EOL MPD Team: NCAR MPD data, Version 1.0, UCAR/NCAR – Earth Observing Laboratory [data set], https://doi.org/10.26023/MX0DZ722M406, 2020. a, b
Nehrir, A. R., Repasky, K. S., Carlsten, J. L., Obland, M. D., and Shaw, J. A.: Water Vapor Profiling Using a Widely Tunable, Amplified DiodeLaserBased Differential Absorption Lidar (DIAL), J. Atmos. Ocean. Tech., 26, 733–745, https://doi.org/10.1175/2008JTECHA1201.1, 2009. a, b, c
Newsom, R. and Sivaraman, C.: Raman Lidar Water Vapor Mixing Ratio and Temperature ValueAdded Products, Tech. rep., DOE Office of Science Atmospheric Radiation Measurement (ARM) Program (United States), https://doi.org/10.2172/1489497, 2018. a, b
NRC: Observing Weather and Climate from the Ground Up, National Academies Press, Washington, D.C., https://doi.org/10.17226/12540, 2009. a
NRC: When Weather Matters, National Academies Press, Washington, D.C., https://doi.org/10.17226/12888, 2010. a
NRC: Weather Services for the Nation, National Academies Press, Washington, D.C., https://doi.org/10.17226/13429, 2012. a
Oh, A. K., Harmany, Z. T., and Willett, R. M.: Logarithmic total variation regularization for crossvalidation in photonlimited imaging, in: 2013 IEEE International Conference on Image Processing, 15–18 September 2013, Melbourne, VIC, Australia, 484–488, https://doi.org/10.1109/ICIP.2013.6738100, 2013. a, b, c, d, e, f, g
Oh, A. K., Harmany, Z. T., and Willett, R. M.: To e or not to e in poisson image reconstruction, in: 2014 IEEE International Conference on Image Processing (ICIP), 27–30 October 2014, Paris, France, 2829–2833, https://doi.org/10.1109/ICIP.2014.7025572, 2014. a
Roelofs, F., Janssen, M., Natarajan, I., et al.: SYMBA: An endtoend VLBI synthetic data generation pipelineSimulating Event Horizon Telescope observations of M 87, Astron. Astrophys., 636, A5, https://doi.org/10.1051/00046361/201936622, 2020. a
Schafer, R. W.: On the frequencydomain properties of SavitzkyGolay filters, in: 2011 Digital Signal Processing and Signal Processing Education Meeting (DSP/SPE), 4–7 January 2011, Sedona, AZ, USA, 54–59, https://doi.org/10.1109/DSPSPE.2011.5739186, 2011. a
Spuler, S. M., Repasky, K. S., Morley, B., Moen, D., Hayman, M., and Nehrir, A. R.: Fielddeployable diodelaserbased differential absorption lidar (DIAL) for profiling water vapor, Atmos. Meas. Tech., 8, 1073–1087, https://doi.org/10.5194/amt810732015, 2015. a, b, c, d, e
Spuler, S. M., Hayman, M., Stillwell, R. A., Carnes, J., Bernatsky, T., and Repasky, K. S.: MicroPulse DIAL (MPD) – a diodelaserbased lidar architecture for quantitative atmospheric profiling, Atmos. Meas. Tech., 14, 4593–4616, https://doi.org/10.5194/amt1445932021, 2021. a, b, c, d, e, f, g, h
Weckwerth, T. M., Wulfmeyer, V., Wakimoto, R. M., Hardesty, M. R., Wilson, J. W., and Banta, R. M.: NCARNOAA LowerTropospheric Water Vapor Workshop, B. Am. Meteorol. Soc., 80, 2339–2357, https://doi.org/10.1175/15200477(1999)080<2331:WOTOFE>2.0.CO;2, 1999. a
Willett, R. and Nowak, R.: Platelets: a multiscale approach for recovering edges and surfaces in photonlimited medical imaging, IEEE T. Med. Imaging, 22, 332–350, https://doi.org/10.1109/TMI.2003.809622, 2003. a
Wulfmeyer, V., Hardesty, R. M., Turner, D. D., Behrendt, A., Cadeddu, M. P., Di Girolamo, P., Schlüssel, P., Van Baelen, J., and Zus, F.: A review of the remote sensing of lower tropospheric thermodynamic profiles and its indispensable role for the understanding and the simulation of water and energy cycles, Rev. Geophys., 53, 819–895, https://doi.org/10.1002/2014RG000476, 2015. a
Xiao, D., Wang, N., Shen, X., Landulfo, E., Zhong, T., and Liu, D.: Development of ZJU HighSpectralResolution Lidar for Aerosol and Cloud: Extinction Retrieval, Remote Sens., 12, 3047, https://doi.org/10.3390/rs12183047, 2020. a, b
Information about the project is available at https://www.arm.gov/research/campaigns/sgp2019mpddemo (last access: 26 August 2022).
 Abstract
 Introduction
 The MPD forward models
 The MPD photoncounting noise model and masking
 The MPD standard method
 The PTVMPD method
 Experiment results
 Conclusion and future work
 Appendix A: Jacobian matrix of MPD forward models (Eq. 5) and gradient matrix of loss function (Eq. 12)
 Appendix B: Poisson negative log likelihood Eq. 8 properties
 Code availability
 Data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Abstract
 Introduction
 The MPD forward models
 The MPD photoncounting noise model and masking
 The MPD standard method
 The PTVMPD method
 Experiment results
 Conclusion and future work
 Appendix A: Jacobian matrix of MPD forward models (Eq. 5) and gradient matrix of loss function (Eq. 12)
 Appendix B: Poisson negative log likelihood Eq. 8 properties
 Code availability
 Data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References