Development of China’s first space-borne aerosol-cloud high-spectral-resolution lidar: retrieval algorithm and airborne demonstration

Abstract

Aerosols and clouds greatly affect the Earth’s radiation budget and global climate. Light detection and ranging (lidar) has been recognized as a promising active remote sensing technique for the vertical observations of aerosols and clouds. China launched its first space-borne aerosol-cloud high-spectral-resolution lidar (ACHSRL) on April 16, 2022, which is capable for high accuracy profiling of aerosols and clouds around the globe. This study presents a retrieval algorithm for aerosol and cloud optical properties from ACHSRL which were compared with the end-to-end Monte-Carlo simulations and validated with the data from an airborne flight with the ACHSRL prototype (A2P) instrument. Using imaging denoising, threshold discrimination, and iterative reconstruction methods, this algorithm was developed for calibration, feature detection, and extinction coefficient (EC) retrievals. The simulation results show that 95.4% of the backscatter coefficient (BSC) have an error less than 12% while 95.4% of EC have an error less than 24%. Cirrus and marine and urban aerosols were identified based on the airborne measurements over different surface types. Then, comparisons were made with U.S. Cloud-Aerosol Lidar with Orthogonal Polarization (CALIOP) profiles, Moderate-resolution Imaging Spectroradiometer (MODIS), and the ground-based sun photometers. High correlations (R > 0.79) were found between BSC (EC) profiles of A2P and CALIOP over forest and town cover, while the correlation coefficients are 0.57 for BSC and 0.58 for EC over ocean cover; the aerosol optical depth retrievals have correlation coefficient of 0.71 with MODIS data and show spatial variations consistent with those from the sun photometers. The algorithm developed for ACHSRL in this study can be directly employed for future space-borne high-spectral-resolution lidar (HSRL) and its data products will also supplement CALIOP data coverage for global observations of aerosol and cloud properties.

Introduction

Aerosols greatly affect the Earth’s radiation budget in both direct and indirectly way through the aerosol-cloud interaction (ACI) [1,2,3]. By analyzing the backscatter laser spectrum or energy variation generated by the interaction between the laser beam and atmospheric particles [4,5,6], lidar has been recognized as a promising remote sensing technique for aerosol and cloud vertical profiling at high spatial and temporal resolution [7,8,9,10]. The ground-based lidar is limited by its spatial coverage, meanwhile the space-borne lidar has great potential to provide global data coverage [9, 11]. China launched the Atmospheric Environment Monitoring Satellite (AEMS) loaded with its first space-borne atmospheric lidar, Aerosol and Carbon Detection Lidar (ACDL) on April 16, 2022 [12, 13] (see Figures S1 and S2). ACDL carries two lidar instruments on a single platform in a 705-km solar synchronous orbit. One is the aerosol-cloud high-spectral-resolution lidar (ACHSRL), and the other is the integrated path differential absorption (IPDA) lidar for atmospheric column CO₂ observations [14,15,16].

The U.S. Cloud-Aerosol Lidar with Orthogonal Polarization (CALIOP) on the Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observations (CALIPSO) satellite, has provided global observations of aerosols and clouds since 2006 [17]. CALIOP has a relatively large uncertainty in extinction coefficient (EC) retrievals of aerosols and clouds as CALIOP was designed as a Mie backscatter lidar system and the lidar ratio (LR) is normally selected or modeled [18]. Taking advantage of the broad spectrum of the Rayleigh scattering from atmospheric molecules, China’s high-spectral-resolution lidar (HSRL) can retrieve the aerosol optical properties and LR directly with a narrow spectral filter to separate the Mie and Rayleigh scattering components from the lidar backscatter signals. The European Earth Clouds, Aerosols and Radiation Explorer (EarthCARE) satellite utilizes a 355 nm Fabry–Pérot interferometer (FPI) HSRL payload, named the ATmospheric LIDar (ATLID) [19], which is scheduled to be launched in 2023. The CALIPSO team has established a data processing system for vertical profiles of aerosols and clouds, including selective iterative boundary locator (SIBYL) for feature detection, scene classification algorithm (SCA), and hybrid extinction retrieval algorithm (HERA) [17]. Meanwhile, the EarthCARE team has developed ATLID level-2 retrieval algorithm (L2a algorithm) using data from EarthCARE simulator. The L2a algorithm mainly covers target detection, aerosol and cloud optical property retrieval and target classification [20, 21].

The ACDL scientific team has studied a variety of aspects about the retrieval algorithm for ACHSRL, starting from ground-based lidar, the error analysis of atmospheric parameters, the optimization design of iodine absorption cell, the calibration of overlap function, and the feature detection to the extinction coefficient retrieval [22,23,24,25,26]. Subsequently, Dong et al. analyzed the optimum iodine absorption line of ACHSRL [27]. Based on Dong’s choice, Liu et al. and Yu et al. have estimated the performance of ACHSRL and showed ACHSRL would be able to provide high quality vertical profiles of aerosols and clouds [13, 28]. Zheng et al. applied the block matching 3D filtering (BM3D) denoising method to ACHSRL attenuated backscatter coefficient (ATB), and the method could meet the noise reduction need [29]. Meanwhile, Mao et al. upgraded the space-borne lidar feature detection using a simple multiscale and minimum cost function method [30, 31].

The airborne ACHSRL prototype (A2P) was developed and an airborne campaign with the A2P was conducted over Qinhuangdao, China, in March 2019, for the evaluation of instrument performance and the validation of retrieval algorithms [32,33,34]. Up to date, the ACDL scientific team is working on the data processing yet for this airborne campaign.

Currently, ACDL is adapting the space environment and will switch on the laser soon. The main purpose of this study is to validate and optimize the corresponding retrieval algorithm of ACHSRL before the satellite observation starts. The Mie signal in the molecular channel of ACHSRL has been eliminated by the iodine absorption cell, resulting in a lower signal intensity compared with the Mie scattering space-borne lidar, which requires special processing in the retrieval algorithm. Hence, the methods of imaging denoising, threshold discrimination and iterative reconstruction are applied in this study to obtain optimal signal-to-noise ratio (SNR), accurate layer distributions and optical properties.

In this study, first, we describe the measurement principle of ACHSRL and its retrieval algorithm. We then present the forward model Monte-Carlo (MC) simulations for a variety of scenes, the airborne campaign in Qinhuangdao, China, and some results from simulations and airborne measurements. The A2P aerosol and cloud optical properties were compared with CALIOP profiles, Moderate-resolution Imaging Spectroradiometer (MODIS) and sun photometers optical depth products. The comparison results showed good agreement between them and robustness of the retrieval algorithm.

Methods

Measurement principle of ACHSRL

Figure 1 is the functional diagram of ACHSRL. The backscatter lights from aerosols and clouds are collected by the telescope and detected by photomultiplier tube (PMT) or avalanche photodiode detector (APD) [35, 36]. A parallel channel and perpendicular channel are designed to detect the linear depolarization of backscatter light while the molecular channel is specialized for the estimate of aerosol and cloud loads.

Fig. 1

Functional block diagram of ACHSRL at 532 nm. ‘F’ = FP etalon filter, ‘L’ = Lens, ‘R’ = Reflector, ‘PBS’ = Polarization Beam Splitter, ‘BS’ = Beam Splitter

The three channels are denoted with subscripts || (parallel), ⊥(perpendicular) and M (molecular). The Mie scattering spectrum is filtered by the narrow-band discriminator iodine absorption cell from the Rayleigh scattering spectrum. Without the strict assumption of the LR, the signals from three channels are sufficient to retrieve the optical properties of aerosols and clouds [13, 37].

The backscatter signal of three channels is expressed as, according to the functional block diagram [13, 37]:

$$P^M(z)=\frac{C^M}{(z-z_0)^2}\lbrack f_m(z)\beta_m^\parallel(z)+f_p(z)\beta_p^\parallel(z)\rbrack T^2(z)$$

(1)

$$P^\parallel(z)=\frac{C^\parallel}{(z-z_0)^2}\lbrack\beta_m^\parallel(z)+\beta_p^\parallel(z)\rbrack T^2(z)$$

(2)

$$P^{ \bot } (z) = \frac{{C^{ \bot } }}{{(z - z_{0} )^{2} }}[\beta_{m}^{ \bot } (z) + \beta_{p}^{ \bot } (z)]T^{2} (z)$$

(3)

$$T^{2} (z) = exp\left[ { - 2\int_{0}^{r} {\alpha (z^{\prime})} dz^{\prime}} \right]{\text{ = exp}}\left[ { - 2\int_{0}^{r} {\alpha_{p} (z^{\prime}) + \alpha_{m} (z^{\prime})} dz^{\prime}} \right]$$

(4)

$$B^{i} \left( z \right) = {{\left[ {P^{i} \left( z \right)\left( {z - z_{0} } \right)^{2} } \right]} \mathord{\left/ {\vphantom {{\left[ {P^{i} \left( z \right)\left( {z - z_{0} } \right)^{2} } \right]} {\left[ {C^{i} } \right]}}} \right. \kern-\nulldelimiterspace} {\left[ {C^{i} } \right]}}$$

(5)

where P is the measured power. C is the channel efficiency constant. The subscripts m and p represent the molecule and particle (aerosols and clouds), respectively. β_m/p(z) is the molecule /particulate backscatter coefficient (BSC), and α_m/p(z) is the molecule/particulate EC at altitude z. T²(z) is the two-way transmittance from the satellite to the scattering volume at altitude z. Bⁱ(z) is the calibrated ATB for three channels (i =||, ⊥, M). f_m is the discriminator transmission of Rayleigh scattering, and f_p is the discriminator transmission of Mie scattering. z-z₀ is the range from the satellite to the scattering volume, z₀ is satellite altitude. The optical properties of aerosols and clouds, such as BSC, EC, LR, depolarization ratio (DR) and optical depth, could be calculated from Eqs. 1, 2, 3, 4 and 5). The derivation is provided in Supplementary Information S1.

Large uncertainty of EC retrieved from this standard method would appear without distinguishing clear air from aerosols or clouds (i.e., feature). Furthermore, negative values of EC could be retrieved due to the low SNR [38]. Therefore, denoising and feature detection methods are utilized and will be illustrated in Retrieval algorithm.

Retrieval algorithm

Figure 2 is the ACHSRL data processing flowchart including the data preparation, pre-processing, and optical properties retrieval. The validation part is illustrated in Results and Discussion.

Fig. 2

The data processing schematic of retrieval algorithm for ACHSRL (A2P)

Data preparation

Data preparation includes the original data, calibration constant and meteorological data, such as temperature, pressure, humidity profiles, GPS, f_m, f_p, β_m, α_m, δ_m, etc. The meteorological data used for ACHSRL retrievals would come from European Centre for Medium-Range Weather Forecasts (ECMWF). During the airborne campaign, most of the meteorological data were obtained from the in-situ instruments on board the aircraft and the nearby radiosondes. Invalid data recorded when the aircraft’s attitude is unstable, i.e. at the takeoff, turning, and landing would be screened out from further analysis. After attitude correction and data screening, the uncalibrated original signal was ready for the pre-processing.

Pre-processing

Pre-processing is aimed at obtaining calibrated ATB. Background noise were first removed by subtracting the mean of the high-altitude clear air signals. Then, the horizontal and vertical averaging were made to guarantee sufficient SNR. According to Eq. 5, the systematic efficiency constant Cⁱ and range factor (z-z₀)² of backscattering signal Pⁱ(z) would be calibrated and turned into ATB Bⁱ(z), which only reflected the backscatter and extinction component. The gain ratio correction of different channels of PMT adopts the Rayleigh signal fitting method to obtain the gain-corrected ATB [39]. In addition, the satellite travels at a speed of 7.5 km/s, resulting in a large difference between adjacent shot signals. Meanwhile, too much average on the original profiles may degrade the true values of the optical parameters. Therefore, BM3D was applied to ATB to preserve layer features since Zheng et al. demonstrated its feasibility in ACHSRL denoising [29]. The output data of pre-processing is the averaged, denoised, and range factor and gain ratio corrected ATB.

Feature detection

Optical properties retrieval in Fig. 2 covers feature detection and EC retrieval. Since the backscatter coefficient ratio (BSR) R(r) is often used to represent the relative strength of aerosols and clouds, BSR was chosen as the threshold determination parameter in the feature detection, benefitting from the direct measurements of backscatter coefficient by HSRL [13]. BSR is defined as:

$$R(r) = 1 + {{\beta_{p} (r)} \mathord{\left/ {\vphantom {{\beta_{p} (r)} {\beta_{m} (r)}}} \right. \kern-\nulldelimiterspace} {\beta_{m} (r)}}$$

(6)

Theoretically, the layers can be detected when R(r) > 1. However, the unavoidable lidar signal noise could lead to fluctuation of R(r) and wrong estimation of vertical layer distribution. Consequently, the SNR should be considered, and the selected R(r) threshold should be dynamically varied with signal noise to prevent the misidentification. The uncertainty of R(r) is relevant with the error of β_p and β_m according to Eq. 6. In this study, the meteorological data for retrieval came from the nearby radiosonde sites, which provided data twice a day. Therefore, the relative errors of f_m, f_p, β_m, α_m and δ_m, which rely on the meteorological data uncertainty, are neglected in the following procedures [26]. Then, the relative error of β_p becomes the critical parameter determining the uncertainty of R(r) and the layer distribution. Mathematical details of the relative error of the backscatter coefficient $\varepsilon_{{\beta_{p} }}$ are shown in Supplementary Information S2 [40]. Assume the dynamical R(r) threshold is T(r), which is expressed as:

$$T\left( r \right){ = }1{ + }{{\varepsilon_{{\beta_{p} }} \left( r \right)} \mathord{\left/ {\vphantom {{\varepsilon_{{\beta_{p} }} \left( r \right)} {\beta_{m} (r)}}} \right. \kern-\nulldelimiterspace} {\beta_{m} (r)}}$$

(7)

the feature would be detected when R(r) > T(r). The threshold method for feature detection could be defined as:

$$\left\{\begin{array}{c}\beta_p\left(r\right)>\varepsilon_{\beta_p}\left(r\right),feature\\\beta_p\left(r\right)<\varepsilon_{\beta_p}\left(r\right),clear\;sky\end{array}\right.$$

(8)

In the following analysis, R(r) was also applied for cloud aerosol discrimination. The R(r) higher than 10 is identified as clouds. Besides, the layers could also be classified as ice clouds with depolarization ratio higher than 0.05, lidar ratio smaller than 40 sr, and temperature lower than -20℃ though the scattering ratio is smaller than 10 [25, 41]. The rest of the layers will be classified as aerosols.

Extinction coefficient retrieval

The iterative reconstruction method was used to retrieve the EC by Xiao et al.[24]. The basic principle of the iterative reconstruction method is a) estimate the close numerical range of LR in different layers based on the Pre-retrieval results from Supplementary Information S1. b) assume the original calibrated ATB and retrieved BSC are accurate. c) repeat the iteration of LR in a different layer to reconstruct the layer ATB. d) compare the reconstructed with original calibrated ATB and select the LR that minimizes the error between the two ATBs. Before the iteration, an approximation of Gaussian probability density function $P_{{S_{P} }}$ was introduced for the lidar signal noise model [42, 43], which is described as:

$$P_{{S_{p} }} { = }\frac{1}{{\sqrt {2\pi } \sigma_{{B^{M} }} }}\exp \left( { - \frac{{\left( {B^{M} - B_{ideal}^{M} \left( {S_{p} } \right)} \right)^{2} }}{{2\sigma_{{B^{M} }}^{2} }}} \right)$$

(9)

For a detected layer, B^M is the molecular channel ATB. S_p is the initial value of the LR for iteration and $B_{ideal}^{M} \left( {S_{p} } \right)$ is the reconstructed ideal molecular channel ATB without noise according to Eq. 5;$\sigma_{{B^{M} }}$ is the standard derivation of B^M. The expected LR would be estimated by minimizing the negative natural logarithm of Eq. 9, which is

$$l\left(S_p,B^M\right)=\left\{\begin{array}{c}\left(B^M-B_{ideal}^M\left(S_p\right)\right)^2/\left(2\sigma_{B^M}^2\right)+{1\mathrm n\left({2\pi\sigma}_{B^M}\right)}^{1/2}\\0\end{array}\begin{array}{c},feature\\,clear\;sky\end{array}\right.$$

(10)

where $l(S_{p} ,B^{M} )$ is assumed to be the loss function. Assume $\overline{S}_{p}$ is the final expected result of the LR, the LR iteration is then turned into a regularization problem:

$$\overline{S}_{p} = \arg \min \{ l(S_{p} ,B^{M} ) + \lambda \cdot ||S_{p} ||_{TV} \}$$

(11)

where arg min{${\mathbf{\mathbb{Z}}}$} means getting S_p that minimizes the ${\mathbb{Z}}$;$\lambda$ is the regularization parameter; The total variation seminorm of the LR $\left\| {S_{p} } \right\|$ is the correspondingly penalty function and constraint for the smoothness of the LR, which could be expressed as [44]:

$$\Vert {S}_{p}\Vert =\sum_{\begin{array}{c}n=1\\ k=1\end{array}}^{\begin{array}{c}N-1\\ K\end{array}}{F}_{n,k}\left|{S}_{p,n,k}-{S}_{p,n+1,k}\right|+\sum_{\begin{array}{c}n=1\\ k=1\end{array}}^{\begin{array}{c}N,\\ K-1\end{array}}{F}_{n,k}\left|{S}_{p,n,k}-{S}_{p,n,k+1}\right|$$

(12)

where, N and K represent that data consists of N range bins on the row axis and K profiles on the column axis; the subscripts n and k are indices to the rows and columns; F_n,k is the layer flag in the row and column (layer would be marked as 1 and clean air would be marked as 0); $\left|{S}_{p,n.k}-{S}_{p,n+1,k}\right|$ and $\left|{S}_{p,n, k}-{S}_{p,n,k+1}\right|$ are the absolute value of LR difference on the row and column direction. Since the signal fluctuation of adjacent profiles in space-borne lidar is larger than that in ground-based lidar, the $\left\| {S_{p} } \right\|$ would be larger compared with ground-based lidar. The regularization parameter $\lambda$ could be smaller for space-borne lidar to constrain the impact of $\left\| {S_{p} } \right\|$ to a reasonable range during the iteration.

Ground-based lidar networks or airborne lidar campaigns are often used for space-borne lidar validation [45, 46]. The original data of ACHSRL would be processed by the satellite research team, which is not available yet. In this study, the forward model simulation of ACHSRL and airborne campaign data in Qinhuangdao, China were used to assess the accuracy and robustness of the retrieval algorithm. The description can be found in Experimental Data.

Experimental data

Forward model simulation

The space-borne lidar forward model generally means the generator of simulated lidar return signals and is often established for performance estimation [47]. The necessary information of ACHSRL forward model includes: the configuration of ACHSRL, provided in Table S1, the noise model for shot noise, solar background noise, dark current noise, thermal noise etc. [42]. The atmospheric scenes are listed in Table 1 [48]. The five individual scenes S1-S5 in Table 1 correspond to low-altitude aerosol, double-layer aerosol, high-altitude aerosol, thin high-altitude cloud, and thick high-altitude cloud. 100 profiles are set in a single scene, and the vertical resolution is 60 m. In our MC simulation, each profile is corrupted by the noise generated randomly based on the noise model. The EC is set to obey exponential distribution at low altitude and Gaussian distribution at high altitude [49]. The averaged extinction coefficient (AEC) is the mean of each layer. DR and LR are constant, and only molecular signals are considered in the clean air.

Table 1 Description of five scenes for the signal simulation of ACHSRL

Based on Eqs. 1, 2, 3, 4 and 5, the MC simulated ATB are shown in Fig. 3. The signal is stronger in parallel channels, whether in the aerosol/cloud layers or the clean air. Meanwhile, the layers are almost filtered out in the molecular channel due to the expected low transmission for the Mie scattering spectrum.

Fig. 3

Simulated ATB of the a parallel b perpendicular, and c molecular channels

Airborne campaign over Qinhuangdao

The A2P was designed to evaluate the performance of the space-borne HSRL and the robustness of the retrieval algorithm. The configuration of A2P at 532 nm is shown in Table S2 [33]. The flight data in Qinhuangdao City, Hebei Province, China, on March 14, 2019 was processed in this study [33, 34, 50]. The A2P flight tracks last from 02:40 to 03:40 (UTC) and cover ocean, coastal city, town, forest, and mountains etc., which is marked in the red curved line in Fig. 4 (a). The green dots in Fig. 4 (a) are the CALIOP laser footprints. For aerosol optical depth (AOD) comparison, two sun photometers (CE318) are facilitated at Beidaihe (119.52°E, 39.85°N) and Funing (119.2°E, 39.9°N) ground station, which are marked in Fig. 4 (a). Also, the Suizhong power plant (120.0°E, 40.1°N) at the coastal area is marked in Fig. 4 (a) with a passing time at 03:14. A2P flew near the Beidaihe and Funing stations at 03:27 and 03:31, and the closest distance was 6.8 km and 14.3 km, respectively. The sun photometers products used for comparison were the AODs at 500 nm wavelength in this study. The original resolution of the A2P signal is 5 m in the horizontal direction (30 Hz repetition frequency, 150 m/s flight speed) and 0.6 m in the vertical direction (250 MHz sampling rate). the averaged A2P signals have a horizontal resolution of 300 m and a vertical resolution of 24 m. To better evaluate the feasibility and robustness of the retrieval algorithm, the flight track is classified into three types based on MODIS surface type products.

Fig. 4

a Map of the A2P flight track in red line and the closest CALIOP laser footprints in green dots. Beidaihe and Funing ground stations with CE318 are marked in asterisks. Power plant is also marker. The background map is the MODIS reflectance in true color observed at 02:25 UTC on March 14, 2019; b Map of MODIS surface types. The flight track over Forest, Town, and Ocean are marked in black, white, and red, respectively; c-d Plots of the CALIOP ATB at 532 nm on Mar 14 and 15, 2019 in color bar

The MODIS surface corrected reflectance is used as a background in Fig. 4 (a) to show the actual atmospheric conditions on March 14 [51]. Based on the MODIS MCD12Q1 surface type product, different surface types on the flight track are displayed in Fig. 4 (b) [52]. The flight track is roughly classified into three surface types: Ocean, Town, and Forest. Town mainly includes urban area, cropland, and grassland in coastal areas; Forest mainly includes farmland, wetland, grassland, shrub, and forest. AOD retrieved by A2P at 532 nm was compared with the MODIS MCD19A2 AOD at 550 nm, which is retrieved with the Multi-Angle Implementation of Atmospheric Correction (MAIAC) algorithm [53]. The overpass time of MODIS was 03:35 (UTC).

CALIOP L2 profiles were also used to compare A2P retrievals [54]. The closest overpass time of CALIOP was 18:25 (UTC) on March 14 and 05:28 (UTC) on March 15, and the closest CALIOP laser footprint is 105 km and 98 km away from the A2P flight track, respectively. The footprints are also grouped into three surface types as shown in Fig. 4 (b). The CALIOP L1 ATB products are shown in Fig. 4 (c)-(d). Cloud layers with the strongest ATB were located at 3–6 km on March 14 and at 10 km on March 15. The profiles of CALIOP located in 38°N -41.5°N were included for comparison.

Results and discussion

Retrieval results with simulation data

BSC, EC, and LR were retrieved for end-to-end simulation and the corresponding errors are shown in Fig. 5. The error of each parameter in Fig. 5 (d)-(f) came from a comparison with the initial assumption in Table 1. From Fig. 5 (a)-(c), the aerosol/cloud layers were accurately identified. The area with largest error appears at 0–3 km in profiles 400–500 mainly due to the rapid laser power attenuation caused by thick clouds at 10–12 km. In the light of the statistical rules for standard deviation, it may be interpreted that 95.4% data of BSC have errors less than 12% and 95.4% data of EC have errors less than 24%. Comparing to the uncertainty of CALIOP official BSC products (~ 20–30%) and EC products (~ 40%), this retrieval algorithm is more consistent and robust [17].

Fig. 5

a-c BSC, EC, and LR retrievals. d-f Error of BSC, EC, and LR retrievals compared with initial assumption in Table 1. The mean and standard deviation of error are shown in the lower left corner

Retrieval results of airborne campaign

The retrieval results of airborne campaign are shown in Fig. 6. The laser path length is estimated from the GPS altitude, the roll and pitch angles of aircraft. Since the laser incident angle varied when the aircraft’s attitude changed, the original signals were removed when the laser path length was much greater than the GPS altitude to avoid invalid retrievals. According to Fig. 4 (b), the flight data have been grouped into seven scene types to evaluate the algorithm performance under different measurement conditions (Ocean1, Town1, Forest1, Town2, Ocean2, Town3, Forest2), separated by red straight lines in Fig. 6.

Fig. 6

Retrievals from A2P measurements on March 14, 2019. a The orange dashed line is the aircraft’s GPS altitude, and the blue line is the laser path length. b-f Depolarization ratio, backscatter coefficient, backscatter coefficient ratio, extinction coefficient, and lidar ratio. g feature classification of clean air, cloud, aerosol and surface

For the A2P surface identification showed in Fig. 6 (g), the ocean altitude is set to be even and close to zero. Since the flight tracks in the Forest area are close to the Yanshan Mountains, the surface altitude in Forest fluctuated up to 0.5–1 km. In the meantime, the Forest surface altitude is close to the CALIOP surface at 40–41°N from Fig. 4 (c)-(d). The above results showed accurate surface detection and range correction in the retrieval algorithm.

There were layers between 2–4 km with BSR higher than 10, LR lower than 30 sr, and drastic scattering change during 02:55 to 03:10 over Forest1 and at 03:34 over Forest2. These layers were detected as clouds and the rest were detected as aerosols [25]. Thick cloud layers were seen at 02:58–02:59, 03:05–03:07, and 03:34 when the laser signals were fully attenuated and no returns from anything below these layers. Meanwhile, semi-transparent clouds were seen at 02:55–02:57, 03:00–03:04, and 03:08–03:10. Note that the DR of clouds at 03:00–03:04 were 0–0.1, lower than others as shown in Fig. 6 (b). These layers could be detected as water or mixed-phase clouds [10, 25]. Over other periods of time, the DR of clouds were higher than 0.15. Since the air temperatures were between -10 and -20 °C in the layers of 2–4 km on Mar 14, 2019, these clouds detected at 02:55–03:00, 03:05–03:10, and 03:34 were probably cirrus (ice) or mixed-phase clouds. The optical properties of cirrus (ice) observed in this study were consistent with those from NASA’s ARCTAS field experiment, which indicated the LR of ice were 10–40 sr and the DR were 0.15–0.4.

Since the selected Ocean area in this study is close to land, the LR of aerosol layers over Ocean1 and Ocean2 ranged from 20 to 40 sr, and their DR were below 0.1 as shown in Fig. 6 (b)-(f). These aerosols may be a mix of marine aerosols with aerosols transported from land, especially along the border of Ocean2 and Town2&3. Such DR and LR of marine aerosols were also previously reported by NASA and DLR’s airborne campaigns (DR from 0 to 0.1 and LR from 10 to 40 sr) [55, 56]. The LR of aerosols over Town and Forest areas were between 30 and 55 sr and their DR varied from 0 to 0.1, suggesting possible urban aerosols and smoke from human activities [55, 56]. As shown in Fig. 6 (c)-(e), the aerosol BSC and EC retrievals over Forest and Town have higher values than those over Ocean. This may be attributed to the emissions from the straw burnings in Forest, the industrial and other human activities in Town. For example, there is a power plant at the coastal area from Fig. 4 (a) [33, 34, 50].

Comparison of A2P profiles with CALIOP

Comparisons with CALIOP L2 aerosol/cloud profiles were made to demonstrate the A2P retrieval accuracy. It may not be appropriate to compare profile-by-profile since the closest distances between CALIOP footprints and A2P tracks are 105 km on March 14 and 98 km on March 15, and the time lag is 15 h and 25 h between two measurements. The BSR of CALIOP was calculated from Eq. 6 and the β_m was estimated from the molecular number density products. The vertical profiles were averaged horizontally for different surface types shown in Fig. 3 (b) for both A2P and CALIOP measurements. The number of averaging are listed in the last column of Table 2. The comparison results between A2P and CALIOP are shown in Fig. 7 and Table 2. Three rows from top to bottom in Fig. 7 are the results over Town, Ocean, and Forest; four columns from left to right are the results of BSR, BSC, EC and LR. Good consistency was found of the profile distributions at 0–3 km over Town, 0–1.5 km over Ocean and 0–4 km over Forest. The BSR, BSC and EC over Forest are much larger than that over Town and Ocean because clouds exist over Forest, while the LRs of A2P distribute more uniformly than that of CALIOP due to the dependence on LR model assumption in CALIOP retrieval algorithm. Statistics of the results are listed in Table 2.

Table 2 Retrieval mean and STD of optical parameters from all the valid bins over Town, Ocean and Forest

Fig. 7

Comparison of BSR, BSC, and EC profiles and statistical probability distribution of LR. a-h The blue (red) solid line is the mean of profiles of A2P (CALIOP) in the along-track direction, and the amplitude of the light blue (red) block is the one order STD. i-l The blue (red) solid line is the mean of aerosol profiles of A2P (CALIOP), while the dashed lines are for clouds. a-d Comparisons in Town; e–h Ocean; i-l Forest

In Table 2, the statistical results are the mean and standard deviation (STD) calculated from all the valid bins over the same surface type. For the comparison over Forest, the aerosol and cloud layers are distinguished to avoid significant bias caused by the data average. The vertical profiles of mean and STD are noted as $\overline{{{\mathbf{\mathbb{Z}}}_{j}^{i} }} \left( z \right)$ and $\sigma_{{{\mathbf{\mathbb{Z}}}_{j}^{i} }} \left( z \right)$, where $\mathbb{Z}=R,\,\beta,\,\alpha,\,S$; i = A2P, CALIOP; j = Town, Ocean, Cloud, Aerosol.

As shown in Table 2, all the STD of the A2P retrievals over Town are smaller than that over Ocean because of the higher SNR from heavy aerosol loads. The LR of A2P over Town is around 35 sr, which suggests the possible mixture of urban and marine in coastal areas. At the same time, the LR over Ocean is smaller (28 sr), probably due to the abundant existence of marine [55, 56]. However, the LR of CALIOP over Town and Ocean are consistent (around 44 sr) mainly attributed to CALIOP’s algorithm. Additionally, the STD of A2P measurements (BSR, BSC, EC) are all smaller than those of CALIOP, especially over Town, indicating the better precision and accuracy of HSRL.

The official L2 vertical feature mask (VFM) products were used for the cloud-aerosol discrimination (CAD) of CALIOP [54], while the A2P clouds were extracted with a BSR > 10 [25]. The profiles were horizontally averaged after the CAD since some clouds and aerosols exist at the same altitude in Forest area as is shown in Fig. 6 (d). The large scattering difference between clouds and aerosols would interfere with the comparison and verification without CAD. The backscatter and extinction intensity of clouds are 1–2 orders of magnitude higher than that of aerosols as shown in Fig. 7 (i)-(k). Due to the spatiotemporal difference, the cloud heights of CALIOP were between 3–4 km, and the layer heights of A2P was between 2.2 and 3.8 km. The STD of A2P cloud heights is lower than that of CALIOP. On the contrary, the STD of A2P aerosol heights is higher than that of CALIOP, which may be influenced by the inaccurate identification of layers (identified as aerosols with 3 < BSR < 10) at the edge of clouds. The cloud LR of A2P is around 22 sr. These high-DR (> 0.15) layers were probably cirrus or mixed-phase clouds, while the aerosol LR of A2P is higher (36 sr) for a mixture of smoke and urban. On the other hand, the LR of CALIOP is much higher than that of A2P, which are 40 sr in clouds and 51 sr in aerosols [55, 56].

The correlations between A2P and CALIOP profiles are shown in Fig. 8. The parameters used for correlation analysis are the Pearson correlation coefficient (R), the mean bias (MB), and the factor of exceedance (FoE). The mathematical definitions are described in Supplementary Information S3 [57]. The correlation coefficients R are 0.79 for BSC and 0.81 for EC in Town, respectively, while the MB are 14.18% for BSC and 21.71% for EC. Moreover, the FoE of BSC in Town is -0.08 (nearly equal), which shows a good agreement with BSC comparison. Since the LR of A2P in Town is lower than CALIOP (34 sr vs 44 sr), most EC of A2P is lower than CALIOP (with a -0.37 FoE). As for comparison in Ocean, the correlation coefficients R are 0.57 for BSC and 0.58 for EC; The MB is 41.13% for BSC and 37.04% for EC. Though the FoE of BSC is zero, discrepancies would still exist due to the lack of available data for CALIOP: only 24 profiles distributed at 0–1.2 km, while the A2P detected 384 profiles at 0–3 km. Except for the data availability, the LR assumption of CALIOP and the spatiotemporal difference between two lidar instruments could also lead to disagreements. In Forest area, the correlation coefficients are 0.82 for BSC and 0.81 for EC. The MB are 140% for BSC and 72% for EC. The aerosol layers significantly affect the large MB at around 3 km from Fig. 8 (f) (h). Besides, the aerosols below the cloud layers often possessed low SNR due to the cloud’s rapid attenuation of laser energy.

Fig. 8

a-b Comparison of $\overline{{\beta_{Town}^{A2P} }} \left( z \right)$&$\overline{{\beta_{Town}^{CALIOP} }} \left( z \right)$, $\overline{{\beta_{Ocean}^{A2P} }} \left( z \right)$&$\overline{{\beta_{Ocean}^{CALIOP} }} \left( z \right)$. The blue (orange) dots are BSC points at different altitudes in Town (Ocean) of A2P and CALIOP; c-d Comparison of $\overline{{\alpha_{Town}^{A2P} }} \left( z \right)$&$\overline{{\alpha_{Town}^{CALIOP} }} \left( z \right)$, $\overline{{\alpha_{Ocean}^{A2P} }} \left( z \right)$&$\overline{{\alpha_{Ocean}^{CALIOP} }} \left( z \right)$. The blue (orange) dots are EC points at different altitudes in Town (Ocean) of A2P and CALIOP; e–f Comparison of $\overline{{\beta_{Cloud}^{A2P} }} \left( z \right)$&$\overline{{\beta_{Cloud}^{CALIOP} }} \left( z \right)$,$\overline{{\beta_{Aerosol}^{A2P} }} \left( z \right)$&$\overline{{\beta_{Aerosol}^{CALIOP} }} \left( z \right)$. The blue (orange) dots are BSC points at different altitudes in Cloud (Aerosol) of A2P and CALIOP; g-h Comparison of $\overline{{\alpha_{Cloud}^{A2P} }} \left( z \right)$&$\overline{{\alpha_{Cloud}^{CALIOP} }} \left( z \right)$, $\overline{{\alpha_{Aerosol}^{A2P} }} \left( z \right)$&$\overline{{\alpha_{Aerosol}^{CALIOP} }} \left( z \right)$. The blue (orange) dots are EC points at different altitudes in Town (Ocean) of A2P and CALIOP

AOD comparison with MODIS and ground-based sun photometers

After cloud screening, the clear-sky A2P AOD retrievals at 532 nm were also compared with MODIS MCD19A2 (550 nm) and sun photometers (500 nm) AOD products, as shown in Fig. 9. The AOD increased from Ocean1 to Town1, Forest1 to Town2, and Ocean2 to Town3, mainly due to the human activities in the Town and coastal area. As expected, the largest AOD change occurred at 03:14 when the aircraft was passing over the power plant, according to Fig. 4 (a). However, the aerosol loads over ocean were relatively low, so there was a decrease of AOD from Town2 to Ocean2.

Fig. 9

Comparisons between A2P AOD at 532 nm and MODIS AOD at 550 nm for cloud-free scenes. a A2P AODs measured from 02:40–03:40 (UTC) are plotted in blue line. MODIS AODs observed at 03:35 (UTC) are plotted in orange line. b Correlation between MODIS and A2P AODs

For better spatial matching between A2P and MODIS data, the mean of the MODIS AOD over the A2P flight track is used for comparison. Consequently, the MODIS AOD in Fig. 9 (a) does not change with time but varies geographically. The comparison in flight time between 03:00–03:30 (UTC) shows that the A2P AODs are overall lower than MODIS AODs. Since the vertical integral path length of A2P (7 km) is shorter than MODIS orbit height (705 km), the MODIS AOD are expected to be higher in cloud-free scenes. The correlation coefficient between the two AODs is 0.71 (Fig. 9 (b)). The FoE is -0.18, which means 68% of A2P AODs are smaller than MODIS AODs.

Table 3 shows the comparison between AODs from A2P and sun photometers (CE318) at Beidaihe and Funing station on March 14, 2019. The time in Column 1 and 5 is the data collection time of sun photometers. The distance between A2P and sun photometers is listed in Column 4 and 8.

Table 3 AOD measurements between A2P and sun photometers (CE318) at Beidaihe (BDH) and Funing (FN) station on March 14, 2019

The diagram of comparison results between A2P and sun photometers is displayed in Fig. 10. The A2P AOD retrievals were generally smaller than those from the sun photometers at both stations except at 02:59 at Beidaihe station. The lower AOD at 532 nm were reasonable after taking into account of the measurement wavelength dependence [58]. The exception at 02:59 may be attributed to the large distance between the A2P and sun photometer. Meanwhile, the mean of AOD at Beidaihe (0.119) is smaller than that at Funing (0.145), which is consistent with the observation from Wang et al. (AODs at Beidaihe and Funing are 0.108 and 0.157) [33]. Overall, the AOD comparison with sun photometer at ground stations showed a good agreement.

Fig. 10

Comparison between AODs from A2P and sun photometers on March 14, 2019. The blue and orange circles are the AODs from the sun photometer at 500 nm band at Beidaihe and Funing station, respectively; the blue and orange squares are the AODs at 532 nm from A2P

Conclusion

This study presented and validated an imperative retrieval algorithm for ACHSRL, which will supplement the U.S. CALIOP data for global lidar observations of aerosols and clouds and will make a significant contribution to atmospheric science research after the successful launch. This study conducted a comprehensive analysis of the algorithm, including ATB calibration, feature detection and EC retrieval using imaging denoising, threshold discrimination, and iterative reconstruction methods.

We demonstrated the algorithm feasibility with the end-to-end MC simulation and showed that the 95.4% of BSC retrievals had errors less than 12% and 95.4% of EC retrievals had errors less than 24%. We also demonstrated the accuracy of feature detection algorithm with the data from the 2019 airborne campaign at Qinhuangdao. Meanwhile, the aerosol classifications were reasonable and consistent. Most aerosols over the Town and Forest regions were identified as mixture from urban sources and smoke from human activities, while over the coastal area or the ocean, the dominant aerosol types are marine aerosols. Furthermore, the cirrus or mixed-phase clouds were also well identified based on their optical properties (BSR > 10, DR > 0.1, LR < 30) and cloud centroid temperature.

Afterwards, the accuracy of the retrievals were validated with satellite observations and ground-based sun photometers. High correlations between A2P and CALIOP vertical profiles of aerosols and clouds were achieved and the errors of A2P BSC and EC were smaller than those of CALIOP except for the aerosol layers over Forest. A2P AODs were consistent with those from MODIS, showing a correlation of 0.71. The spatial variation of A2P AODs results was reasonable over a variety of surface types.

In the future, the A2P tracks could fly along the footprint of ACHSRL to obtain a better spatiotemporal correlation and allow better comparison and validation. The accuracy of the algorithm can be further validated and optimized for variety of atmosphere scenes and surface conditions with the measured data from ACHSRL in the near future.