Two-dimensional and multi-channel feature detection algorithm for the CALIPSO lidar measurements

In this paper we describe a new two-dimensional and multi-channel feature detection algorithm (2D-McDA) and demonstrate its application to lidar backscatter measurements from the Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observations (CALIPSO) mission. Unlike previous layer detection schemes, this context sensitive feature finder algorithm is applied to a 2D lidar "scene"; i.e., to the image formed by many successive lidar profiles. Features are identified when an extended and contiguous 2D region of enhanced backscatter signal rises significantly above the expected "clear air" value. 5 Using an iterated 2D feature detection algorithm dramatically improves the fine details of feature shapes and can accurately identify previously undetected layers (e.g., subvisible cirrus) that are very thin vertically but horizontally persistent. Because the algorithm looks for consistent 2D patterns, it potentially offers improved discrimination of juxtaposed cloud and aerosol layers. Moreover, the 2D detection algorithm uses the backscatter signals from all available channels: 532 nm parallel, 532 nm perpendicular, and 1064 nm total. Since the backscatter from some aerosol or cloud particle types can be more pronounced 10 in one channel than another, simultaneously assessing the signals from all channels greatly improves the layer detection. For example, ice particles in subvisible cirrus strongly depolarize the lidar signal and, consequently, are easier to detect in the 532 nm perpendicular channel. Use of the 1064 nm channel greatly improves the detection of dense smoke layers, because smoke extinction at 532 nm is much larger than at 1064 nm, and hence the range-dependent reduction in lidar signals due to attenuation occurs much faster at 532 nm than at 1064 nm. Moreover, the photomultiplier tubes used at 532 nm are known to 15 generate artifacts in an extended area below highly reflective liquid clouds, introducing false detections that artificially lower the apparent cloud base altitude, i.e. the cloud base when the cloud is transparent or the level of complete attenuation of the lidar signal when it is opaque. By adding the information available in the 1064 nm channel, this new algorithm can better identify the true apparent cloud base altitudes of such clouds. 1 https://doi.org/10.5194/amt-2020-369 Preprint. Discussion started: 29 September 2020 c © Author(s) 2020. CC BY 4.0 License.


Introduction
The Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observation (CALIPSO) mission (Winker et al., 2010) has provided direct measurements of cloud and aerosol vertical distributions with a very high vertical resolution since 2006. A key component of these measurements is made by the active remote sensing instrument CALIOP (i.e., the Cloud-Aerosol Lidar with Orthogonal Polarization), a two-wavelength (532 nm and 1064 nm), polarization sensitive elastic backscatter lidar. 25 The knowledge of the cloud and aerosol vertical distributions and their properties is critical in assessing the planet's radiation budget (e.g., Shonk and Hogan, 2010), in evaluating the atmospheric radiative heating rate (e.g., Huang et al., 2009), and for advancing our understanding of cloud-climate feedback cycles that occur as the climate warms (e.g., Tsushima et al., 2006).
The critically important first step in retrieving the spatial and optical properties of clouds and aerosols is to determine where these "features" are located in the vertical, curtain-like images (altitude vs. satellite track) of the backscattered lidar signals the contextual information provided by a continuous series of profile measurements by searching for cloud and aerosol patterns in the two-dimensional (2D) image formed by successive lidar profiles. If edge detection techniques based on 2D gradient 45 search routines do not seem adequate for spatial analysis of lidar data (Vaughan et al., 2005), methods based on sliding window operations have shown to greatly improves the feature shape detection (e.g., Hagihara et al., 2010;van Zadelhoff et al., 2011;Herzfeld et al., 2014).
Here we propose a new 2D and multi-channel feature detection algorithm (2D-McDA). This "context sensitive" feature finder algorithm is then applied to a 2D lidar signal scene; i.e., to the image formed by many successive lidar profiles. Moreover, the 50 2D detection algorithm uses the backscatter signals from all available channels: the 532 nm co-polarized (or parallel) signal, the 532 nm cross-polarized (or perpendicular) signal, and the 1064 nm signal. Since the backscatter from some aerosol or cloud particle types can be more pronounced in one channel than another, simultaneously assessing the signals from all channels is expected to greatly improve the layer detection.
Section 2 presents a refined method for determining feature detection thresholds, which are a critically important component 55 of the detection algorithm. Section 3 presents the detection algorithm. The detection of the Earth's surface is described first as it is performed first and separately from the cloud and aerosol detection. This has been shown to have many practical advantages.
Then, the cloud and aerosol detection algorithm is described. Finally, the detections from each channel are merged into a composite feature detection mask. Section 4 shows how this new algorithm improves the feature detection compared to the CALIPSO version 4 vertical feature mask (VFM). 60 2 Threshold based feature detection Atmospheric lidars measure attenuated signal backscattered by molecules (m) and particles (p) β (r) = (β m (r) + β p (r)) T m (r) 2 T O3 (r) 2 T p (r) 2 , where β m (r) and β p (r) are the volume backscatter coefficients for molecules and particulates, and T m (r) 2 , T O3 (r) 2 and T p (r) 2 are, respectively, the two-way transmittances for molecules, ozone, and particles, and r is the range from the satellite altitude. If there are no particles in the atmosphere, Eq. (1) reduces to the molecular attenuated backscatter coefficient A feature, i.e., a cloud or an aerosol layer, appears as an extended and contiguous region of enhanced attenuated backscatter signal that rises significantly above the expected clear sky (molecules only) value. To detect these features, a first step is to determine a threshold above which the attenuated backscatter signal should be considered as due to an atmospheric feature. Not 70 all signals larger than the expected molecular attenuated backscatter signal β m (r) are necessarily due to real features; instead they can be caused by noise. A range-dependent threshold is computed following where ∆β b (r) is the standard deviation due to background noise (range-independent 1 ), ∆β m (r) is the expected standard deviation due to the shot noise (range-dependent) in the expected clear sky, and N (r) is the number of bins averaged. The 75 f corr (r) term is a correction function which takes into account the partial vertical correlation in samples due to the limited electronic bandwidth and the shifting and rebinning that can occur in the altitude registration phase of the level 1A processing (details in Appendix A). The number of shot noise standard deviations considered in the threshold is quantified by the factor k, which can be tuned according to the degree of sensitivity needed to avoid false detections. β m (r) is derived from modeled profiles of molecular and ozone number densities. ∆β b (r) is derived from the on-board computation of the root mean square 80 (RMS) of the background signal in the high altitude background region (HABR) between 65 and 80 km for each shot . ∆β m (r) is estimated using its proportional relation with the square root of β m (r) (e.g., Liu and Sugimoto, 2002), The NSF is evaluated from the solar background signal during daytime for the 532 nm parallel and perpendicular channels 85 Liu et al., 2006). At 1064 nm, CALIOP uses an avalanche photodiode (APD) detector, rather than the photomultiplier tubes (PMTs) that are used for the 532 nm channels. Because the APD dark noise overwhelms the 1064 nm shot noise, only the background noise is considered at 1064 nm. Figure 2 shows the range-dependent threshold (red) computed from Eq. (3) with k = 2 applied to the 532 nm parallel lidar signal (blue) for a clear-sky case study during nighttime. Note the noise due to the quantum nature of photons (shot-noise) in 90 this figure. Indeed, although background noise, mainly due to solar radiation, is quite low during nighttime, the lidar signal shows large variations around the expected clear-sky return (black). The range-dependent threshold correctly keeps most the signal below the detection level. Jumps at 8.2 and 20.2 km reveal the change of onboard averaging resolution. However, a few points of the lidar signal still exceed the threshold. Some continuity tests are then needed to determine whether these high signals are due to noise or instead part of an extended feature. Unlike the current CALIPSO detection algorithm, this continuity 95 test will be applied in two dimensions. Figure 3 shows all pixels of Fig. 1 where the lidar signal is above the range-dependent threshold computed from Eq. (3) with k = 2.
For practical reasons, the layer detection scheme uses the attenuated scattering ratio rather than the attenuated backscatter. The attenuated scattering ratio threshold is then obtained from Equation (6) is applied to the three lidar channels (532 nm parallel, 532 nm perpendicular, and 1064 nm) during the 2D-McDA process.

2D and multi-channel feature detection algorithm
2D-McDA is applied to the scattering ratio signals at 532 nm parallel, 532 nm perpendicular, and 1064 nm. First, the detection 105 of the surface altitude is performed and the signal from this altitude and below is removed from the data. Second, the detection of cloud and aerosol layers is done in each channel based on iterated detection thresholds and 2D spatial continuity tests.
Finally, the masks from the three channels are merged in a composite feature mask.

Surface detection
Before detecting cloud and aerosol layers using the detection threshold as described in the previous section, we perform first 110 an independent surface detection. Doing the surface detection in a first and separate step allows a better retrieval of the surface echo and prevents complications in the cloud and aerosol layer detection process. Also, knowing where the surface is detected allows subsequent separation of semi-transparent features from opaque features, which is essential for accurately estimating range-resolved profiles of extinction coefficients (Young et al., 2018). Operationally, atmospheric features are defined as being opaque when no surface return or other atmospheric feature can be detected below them. From this definition it follows that the 115 signals received from beneath opaque features have been fully attenuated within these features. The surface detection algorithm used here is closely akin to the one described in Vaughan et al. (in progress) and is applied to the 532 nm parallel and 1064 nm channels (details in Appendix B). The signals from the top of the detected surface echo and below this point are removed from the data. To minimize computation times, the surface detection algorithm is not applied to the 532 nm perpendicular channel signal. The backscatter from ocean surfaces (covering~70% of the planet) does not depolarize and, excluding snow and ice, the 120 depolarization of most land surfaces is relatively low (Lu et al., 2017), hence the preponderance of the surface backscatter is in the parallel channel. The altitude retrieved in the parallel channel is used to remove signal at and below the estimated surface altitude in the perpendicular channel. Note that there is some small chance that a surface echo can appear in the perpendicular channel but not be visible/detected in the parallel channel. The detection of the surface corresponding to Fig. 1 is shown in brown on Fig. 3.

Cloud and aerosol detection
The detection of cloud and aerosol layers in a single channel curtain of lidar measurements takes place in four main steps: 1. Detecting strong features; i.e., identifying contiguous regions of enhanced attenuated scattering ratios that rise above the feature detection threshold (which is repeatedly decreased from very large value of k down to k = 1) without applying any signal averaging (i.e., d = 1 − 4 in Table 1, refer to Sect. 3.2.1); 130 2. Flagging regions below opaque features as "fully attenuated" (FA) and regions below transparent features where the signal is strongly attenuated with the low confidence flag "almost fully attenuated" (AFA) (Sect. 3.2.2); 3. Averaging of those signals not already flagged using a horizontal sliding window (Sect. 3.2.3); 4. Detecting faint features; features are once again identified as contiguous regions of enhanced signal (i.e., averaged attenuated scattering ratios) that rise above the recomputed feature detection threshold (Sect. 3.2.1).
135 Figure 4 shows the flowchart of the whole detection algorithm. The parameter values used at the different detection levels are given in Table 1.
The following subsections give the details of the main steps presented above.

Load level 1 profile data
Surface detection

Detection
The detection phase is performed following three substeps: 140 1. All pixels of within image that exceed the threshold are first flagged as detected (Fig. 3); 2. A spatial coherence test window is applied to the image of detected/undetected pixels. It smooths the shape of detected pattern and remove isolated noisy detected pixels by turning some of detected pixels to undetected or undetected to detected; 3. Smoothed patterns are required to meet a minimum numeric threshold of contiguous pixels. Patterns that fail to meet this 145 threshold are removed from consideration for this level of detection.
The scattering ratio image used in the layer detection scheme has a spatial resolution of one laser pulse horizontally and 30 m vertically, equivalent to the finest spatial resolution of the CALIOP data. As described in Hunt et al. (2009), CALIOP data is averaged onboard the satellite with spatial resolutions that vary according to altitude. Scattering ratios in regions where the data resolution is coarser than the image resolution are duplicated to the image resolution. The first substep is then to flag all 150 pixels of this image which exceed the detection threshold given by Eq. (6) with the value of k defined in Table 1. For example, for the 532 nm parallel channel, at the detection level d = 3, the detected pixels are those where signal is greater than 1 + k times the expected noise standard deviation with k = 2 ( Fig. 5a; orange pixels). Then, the second substep is to apply a spatial coherence test window (rows a in Table 1) on these detected pixels ( Fig. 5a,b). Here, an 11×11 pixels window is applied to each pixel of the image, with the window being centered successively on all pixels. If the number of originally detected pixels 155 in the window is greater than half of the total number of pixel in the window (≥ 61 for a 11×11 pixels window), then the center pixel is considered as detected. If not, the center pixel is considered as not detected. In this smoothing step, the determination of detection status does not rely on a single pixel exceeding its threshold, but instead on the fraction of neighboring pixels that exceed their thresholds. Consequently, a pixel classified as detected may not itself exceed the detection threshold. Similarly, a pixel that exceeds the threshold may not ultimately be classified as detected. The pixel count within the window is limited 160 to those detected at the current detection level d and at the previous detection level d − 1. This allows detection continuity of  Other flagged pixels (i.e., "Surface", detection ≤ d − 2, "Likely artifact", "Fully attenuated", "Almost fully attenuated", and "Low confidence small strips" (to be described in detail later in this subsection and subsection 3.2.2)) are not considered in the window and the total number of candidate pixels in the window is decreased accordingly. The shapes of detected features 165 are smoothed by the spatial coherence test window, while the noise (isolated orange pixel) is removed. However, some small clusters of pixels sometimes persist. Those small clusters cannot be confidently declared as feature at this stage. They can be due to noise or they can be part of a larger, fainter feature. Then, we decide not to consider these small patterns as detected feature and retain these regions for inclusion in the signal averaging used in successive iterations of the algorithm. To be declared as detected feature, smoothed detected patterns need to consist of more than n connected pixels (see Table 1), otherwise they are 170 removed from this detection level.
This detection procedure is applied several times with different thresholds, different spatial coherence test window, and different limits on the number of connected pixels required (Table 1) in order to detect all layers from the most evident, very strong patterns to the very faint ones, and from geometrically small patterns to very extended ones. Note that the horizontal spatial coherence test window (3×21) enables the detection of faint but horizontally extended cirrus such as the layer shown 175 between 50 km and 100 km in Fig. 1. The detection of this subvisible cirrus is presented in Fig. 6. Figure 6a to Fig. 6b shows the implementation of the 3×21 spatial coherence test window. We see that the cirrus pattern is smoothed and now clearly 9 https://doi.org/10.5194/amt-2020-369 Preprint. Discussion started: 29 September 2020 c Author(s) 2020. CC BY 4.0 License. appears on Fig. 6b due to the fact that most of the noise around has been removed. However, many small clusters of pixels persist. By applying the minimum numeric threshold on the detected pattern, we are able to remove small cluster due to noise while keeping the real cirrus (Fig. 6c).

Special flags
For the 532 nm channels, a first detection of very strong signal is performed (see d = 1 in Table 1). The aim of this initial scan is to identify the tops of very strongly scattering liquid clouds and ice clouds containing high fractions of horizontally oriented ice (HOI) crystals. The non-ideal transient response by PMTs following these very strong signals often generates a spurious, exponentially decaying signal enhancement in the underlying range bins (McGill et al., 2007;Hunt et al., 2009;Lu et al., 2020).

185
The presence of these "noise tails" in the 532 nm signals can introduce large biases into the determination of the apparent bases of opaque water clouds. To exclude this artifact in the detection process, the 600 m below the base of the detected very strong signal are flagged as "Likely artifact" and removed from the signal. Since the APD used in 1064 nm channel does not produce these noise tails, we rely on the 1064 nm channel for the detection of the apparent base of these strongly scattering layers.
After detection of the strongest features, i.e without signal averaging (d = 1 − 4 in Table 1), we flag as "Fully attenuated" After removing all data identified with these low confidence flags from the attenuated scattering ratios, the signal is averaged in order to try to detect fainter features.

Signal averaging
We then average the remaining signal (here the attenuated scattering ratios) using an edge-preserving Gaussian sliding window that extends over 5 km (15 profiles) horizontally and a single range bin vertically (a in Table 1). A Gaussian weight with a standard deviation of 1.67 km is applied, thus giving a stronger weight to pixels closer to center of the window than at the edges. We chose a horizontal window here because the spatial extent of very faint layers is mainly in the horizontal direction.

205
Once the averaging is performed, the detection substeps (Sect. 3.2.1) are then applied to the averaged signal. Note that any feature split by a low confidence vertical band (FA, AFA, and small strips) is considered as a unique feature when counting the number of connected pixels. Figure 7 shows the final mask for the 532 nm parallel channel after the detection of the faint features.

Three channels composite detection 210
The detection algorithm is applied individually to the lidar signal from each of the three channels (Fig. 8) and all pixels identified as features in any of the three channels are retained in the composite mask (Fig. 9a). Comparing this new feature mask (Fig. 9a,b) to the current version of the Vertical Feature Mask (VFM) (Fig. 9c), we first note the improvement in the detected contour of the large cirrus. We also note that the 2D-McDA readily detects faint cirrus (e.g., as seen between 0 km and 75 km) that is missed by the current VFM. The vertical spreading of the clouds seen in the VFM at around 7.5 km in altitude 215 and between 500 km and 900 km horizontally is due to the aforementioned PMT artifact afflicting the 532 nm signals beneath strongly scattering layers. This is not seen in the 2D-McDA feature mask because pixels below the cloud top are flagged as "Likely artifact" in the 532 nm channels and so we make no attempt to retrieve the cloud apparent base of such opaque clouds at this wavelength. Instead, in these cases we retrieve the true penetration depth estimates using the 1064 nm signals (Fig. 8c), which are not affected by detector transient response artifacts (see light blue = "1064 only").

4 Performance assessments and comparisons to version 4
In this section, we present two case studies to show the improvements bring by this new feature detection approach. Figure 10 presents a dense smoke event in Siberia on July 26, 2006 during daytime. The smoke layer is opaque at 532 nm and thus we do not see any surface echo for this channel (Fig. 10a). Note that the smoke is non-depolarizing so there is no 225 perpendicular signal (Fig. 10b). Because the standard CALIOP layer detection only examines the 532 nm channel, the VFM ( Fig. 11c) indicates that the signals are fully attenuated after detecting (at 532 nm) the apparent base of the smoke layer.

Dense smoke
However, at 1064 nm the dense smoke layer is transparent and the surface is readily detected (Fig. 10c). This scene clearly illustrates the advantage gained by using a multi-channel feature detection algorithm, since the full vertical extent of the smoke plume can only be retrieved by inspecting the 1064 nm measurements (light blue in Fig. 11a).  Figure 12 presents the attenuated backscattered lidar signal in the three channels for another case study showing a variety of cloud types and shapes which occurred above Ethiopia on August 31, 2018 during nighttime. We can see that the artifacts below liquid water clouds (close to the surface and up to 8 km) appear in the 532 nm parallel (Fig. 12a) and the 532 nm perpendicular ( Fig. 12b) channels but not at 1064 nm (Fig. 12c). We note that thin cirrus clouds, like the one at 17 km in altitude between 235 1550 km and 1850 km, are clearly brought out in the 532 nm perpendicular channel (Fig. 12b). If we look now at the composite feature detections derived from these three signals (Fig. 13a), we note again how well the apparent bases of liquid clouds are retrieved by using the 1064 nm channel. We note also that the successful identification of thin cirrus can largely be attributed to our use of the 532 nm perpendicular channel. Figure 13b shows the same mask as Fig. 13a but with the same colors that are used for the VFM images (Fig. 13c). This change of colors is intended to facilitate one-to-one comparisons between the two detection schemes. However, note that the yellow and white colors do not discriminate aerosol from cloud, as in the VFM, but instead simply differentiate weak from strong features based on whether the feature detection required data averaging (yellow) or not (white). Finally, Figure 13d shows the difference between the new composite feature detection mask and the VFM.

Variety of cloud type and shape
We see that the contour of features retrieved by the 2D-McDA represents a distinct improvement over the squared boundaries reported by the VFM. We note too that the new algorithm detects thin clouds that are obviously missed by the VFM and that it 245 eliminates significant detection artifacts reported by the VFM between 700 km and 900 km.

Conclusions
This

255
-Because it is applied to single profiles averaged over several different horizontal resolutions, the standard CALIOP feature detection produces blocky, rectangular layers. The complex shapes of aerosol and cloud features are better preserved by the 2D-McDA windowing and data aggregation operations, which provide the flexibility required to distinguish fine spatial details. It is hoped that this improved feature detection will lead to improvements in classifying features according to type (e.g., clouds vs. aerosols) and in their optical property retrievals. Ideally, it will also offer improved discrimination 260 between juxtaposed cloud and aerosol layers or identifiable regions of ice and liquid water in a cloud. The improvement of the cloud shape detection is by itself important for example for studies interested in anvil clouds (e.g., Bony et al., 2016;Hartmann, 2016); -The detection of subvisible cirrus is significantly enhanced by both the 2D detection scheme and the use of the 532 nm perpendicular channel, which is especially sensitive to the presence of depolarizing ice crystals. Those clouds play an im-265 portant role in the climate system as they regulate the vertical transport of water vapor near the upper troposphere-lower stratosphere (e.g., Jensen et al., 1996;Luo et al., 2003), influence the local thermal budget, and drive dynamics of the tropopause region (e.g., Hartmann et al., 2001;McFarquhar et al., 2000); -The apparent base altitudes of highly reflective clouds, i.e. the levels of complete attenuation of the lidar signal, which are routinely biased low (by several hundred meters) due to the non-ideal transient response of 532 nm photomultiplier tubes, smokes, the 1064 nm signals are attenuated significantly less than at 532 nm and hence can generally penetrate the full vertical extent of these layers. Those biomass burning aerosols play a significant role in the Earth's radiative balance by their scattering and absorption of incoming solar radiation (e.g., Penner et al., 1992;Christopher et al., 1996) and the interaction they have with clouds (e.g., Kaufman and Fraser, 1997). The full detection of those layers will lead to more accurate aerosol optical depth retrievals, which will improve estimates of the radiation budget and interactions with where Cov(·, ·) is the covariance. When averaging N bin consecutive 15-m bins, we can consider they have approximately the same range and then that their variance is constant: Var(b i ) = Var(b). If the b i bins were uncorrelated, then we would have 295 13 https://doi.org/10.5194/amt-2020-369 Preprint. Discussion started: 29 September 2020 c Author(s) 2020. CC BY 4.0 License.
Cov(b i , b j ) = 0, ∀(i = j), and then Var(B) = Var(b) N bin . However, since each bin is partially correlated with its vertical neighbors, we have Cov(b i , b i+m ) = constant ∀i for each lag of m range bins. Then, Eq. (A1) can be rewritten following where R(m) = Cov(bi,bi+m) is the autocorrelation coefficient for a lag of m range bins. It follows that the correction function to apply on the standard deviation to take into account the vertical partial correlation due to the electronic bandwidth is 300 f a (N bin ) = 1 + 2

A2 Correction due to redistribution in altitude registration
The altitude for each sample of a raw data profile is recalculated with more accurate information about the satellite altitude and laser viewing angle in the data processing on ground. In general, shift for a few range bins is needed for the full resolution (30 m) samples. Redistribution and rebinning of two neighboring samples whose range resolutions are coarse (60 m, 180 m, 305 and 300 m) are performed, which can induce additional correlation to the redistributed samples that also needs to be taken into account when determining the detection threshold. When there is a shift of N shift15 15-m bins, each new shifted bin B k is computed from the weighted average of the two original bins (with N bin size resolution) it steps across B k and B k+1 (Fig. A1) following and It follows that the correction function to apply on the standard deviation to take into account both the vertical partial correlation due to the electronic bandwidth and the redistribution in altitude registration is Appendix B: Surface detection 320 The aim of this procedure is to detect a surface echo in the near neighborhood region of the estimated surface altitude z surf given by a digital elevation model (DEM). The width of this region will vary according to surface type. Since we are highly confident of the surface altitude over the ocean (where z surf = 0), we will only search in a very narrow range of altitudes for profiles measured over the ocean. On the other hand, we are somewhat less confident of the DEM surface altitudes reported over land and even much less confident over Greenland and Antarctica, so our search regions over land will be larger. The 325 surface echo can be very weak due to attenuation by aerosol and cloud layers above. Then, we try to detect even the weakest surface echo as long as it is substantially above background noise. This procedure is applied at single-shot resolution only. For each shot, the method is made up of the following steps: 1. Compute r surf , the estimated range of the surface, from z surf and the satellite altitude z sat ; 2. Compute ∆β b ( r surf ), the standard deviation due to background noise in the β domain at the range r surf ; (a) If surface type is Water and z surf = 0, then ∆i = 2 (≡ 60 m), 6. Determine z min and z max , the altitudes of the minimum and maximum values of the derivatives in the surface search region, and i min and i max , their respective bin index; 340 7. Determine β max , the maximum signal magnitude lying between z min and z max 2 ; 8. Sequentially test the three following rules to determine if we have identified a legitimate surface return: (a) z min > z max ,  . Before these substeps, surface is first detected (brown), then very strong signal (k = 100) occurring on highly reflecting liquid clouds is detected (black) and the 600 m below is flagged as "likely artifact" (gray) as it is the region where we see artifacts due to the time response of photomultiplier tubes (PMTs) in the 532 nm channels. Two detections were also made: one with k = 20 and another with k = 2, a 11×11 spatial coherence test window, and n = 60.