A Coupled BRDF CO 2 Retrieval Method for the GF-5 GMI and Improvements in the Correction of Atmospheric Scattering

: The Greenhouse Gases Monitoring Instrument (GMI), on board the Chinese Gaofen-5 (GF-5) satellite, provides rich observation data for the global remote sensing of atmospheric CO 2 . To meet the high-precision satellite retrieval needs of atmospheric CO 2 , this paper designs a coupled bidirectional reﬂectance distribution function (BRDF) CO 2 retrieval (CBCR) method, which describes the surface reﬂectance characteristics by the BRDF, corrects for atmospheric scattering based on full physics retrieval theory, and ensures the stable retrieval of multiple parameters and atmospheric CO 2 by enriching prior constraints. Theoretical analysis shows that the inﬂuence of atmospheric scattering induced by the surface bidirectional reﬂectance characteristics is signiﬁcantly related to the aerosol optical depth (AOD), solar zenith angle (SZA), and viewing zenith angle (VZA). The validation of GMI CO 2 retrievals shows that the CBCR method signiﬁcantly reduced the inﬂuence of the surface bidirectional reﬂectance characteristics under high AOD and high SZA conditions, decreased the atmospheric CO 2 retrieval error from 0.58 ± 5.64 ppm to − 1.33 ± 3.13 ppm, and increased the correlation with the temporal variation of actual atmospheric CO 2 from 34.7 to 76.8%. Our CBCR method can correct the inﬂuence of atmospheric scattering induced by the surface bidirectional reﬂectance characteristics on atmospheric CO 2 retrieval, and this work demonstrates that describing the surface reﬂectance characteristics by using BRDF is a promising idea in the ﬁeld of satellite CO 2 retrievals.


Introduction
Atmospheric carbon dioxide (CO 2 ) is the most important anthropogenic greenhouse gas. The globally averaged CO 2 concentration has rapidly increased over the last 250 years and recently reached approximately 400 ppm; this increasing concentration has greatly contributed to the global radiative force [1,2]. Hence, global measurements of the atmospheric CO 2 concentration are remarkably important for better understanding the Earth's carbon cycle and CO 2 sources and sinks and for assessing future trends associated with global warming.
Remote sensing with satellite sensors is one of the most effective approaches for providing dense and uniform atmospheric CO 2 observations on a global scale, due to their global coverage and high sampling frequency, which can better meet the application needs of climate research. Accordingly, satellite remote sensing technology for atmospheric CO 2 has developed rapidly in recent years. To date, several satellites have been deployed to monitor greenhouse gases, such as the Greenhouse Gases Observing Satellite (GOSAT) series, the Orbiting Carbon Observatory (OCO-2) series, Sentinel-5P, TanSat, FengYun-3D (FY-3D), and GaoFen-5 (GF-5) [3].
The satellite retrieval accuracy is critical for application purposes. Satellite measurements with the accuracy of approximately 1% (~4 ppm) for 1000 × 1000 km 2 regional averages and monthly means have the potential to significantly enhance the accuracy of CO 2 source and sink estimations [4,5]. However, the environment of satellite observations is complex, and their retrieval accuracy suffers from the combined effects of the atmosphere, surface, and instruments used. Several studies have shown that the largest error source is atmospheric scattering. Consequently, to achieve a high retrieval accuracy, many correction schemes have been designed for CO 2 satellite remote sensing retrievals, particularly for the atmospheric scattering corrections of GOSAT and OCO-2, and these schemes have been optimized to obtain a greater number of effective products and enhance the ability to quantify carbon sources and sinks [6][7][8].
The Chinese GF-5 satellite, which can provide an abundance of data for the global detection of atmospheric CO 2 , was successfully launched at the Taiyuan Satellite Launch Center on 9 May 2018. The Greenhouse Gases Monitoring Instrument (GMI) equipped on GF-5 was designed to monitor the global distributions of greenhouse gases (namely, CO 2 and CH 4 ) from space. After an initial period of heating to remove the effects of H 2 O contamination, GMI began to work on 29 May 2018 [9]. GF-5 travels within a sun-synchronous orbit at an altitude of 705 km, with an equatorial crossing local time descending node (LTDN) of 13:30. This satellite completes an orbit in approximately 100 min and operates on a global basis with a 7 day orbit repeat cycle but with a 51-day footprint revisit cycle, and provides data with a spatial resolution of 10.5 km. The common measurement mode of GMI is the nadir five-point cross-track scan mode. The footprints are separated by approximately 230 km in the cross-track direction and 102 km in the along-track direction at 30 • latitude in the northern hemisphere, which is similar to GOSAT.
To retrieve high-precision CO 2 products from GMI measurements, we improved the initial retrieval algorithm, named GMI_XCO 2 , which is based on the principles of differential optical absorption spectroscopy (DOAS), does not correct for atmospheric scattering, and uses directional polarimetric camera (DPC) aerosol products to screen data with an aerosol optical depth (AOD) of less than 0.3 for the retrieval of CO 2 contents [10]. Theoretically, the initial method achieves a good retrieval accuracy, but when it is applied to actual retrieved GMI observation data, we found that the retrieval results are highly discrete, and few data can reach an accuracy of 1%. According to research on CO 2 remote sensing retrievals from the measured data of GOSAT and OCO-2, the influence of atmospheric scattering on the high-precision CO 2 retrievals via satellite remote sensing is complex and important, and therefore, it must be corrected [11,12]. Likewise, the bidirectional reflectance distribution function (BRDF) characteristics of the surface [13][14][15] and the corresponding sensitivity to atmospheric scattering cannot be ignored [16,17]. Therefore, based on the theory of the full physics (FP) algorithm [18,19], we designed a conventional FP method to synchronously invert the surface reflectance, aerosols and CO 2 to correct the influence of atmospheric scattering, and on this basis, we further considered the influence of surface BRDF characteristics to realize the high-precision retrieval of GMI-measured data. We call this method the coupled BRDF CO 2 retrieval (CBCR) method, and we applied it to the retrieval of GMI-measured data in this study.
In this paper, we present the designed CBCR method, which improves the satellite retrieval ability as a result of the coupling of the BRDF model with a user-defined aerosol model and the enhanced control of a priori constraints. The retrieval of GMI measurements above Total Carbon Column Observing Network (TCCON) stations was carried out, the retrieval results of the CBCR method are displayed, and the improvement in the retrieval accuracy over the original algorithm was analyzed. The remainder of this article is organized as follows: Section 2 briefly introduces the GMI equipment and data. Section 3 introduces the CBCR method. In Section 4, the influence of the surface bidirectional reflectance on the retrieval accuracy is analyzed from a theoretical perspective, the two methods are compared, and the main influencing factors are delineated. Section 5 compares the two retrieval methods using GMI measurements and analyzes the correction effect of the new Remote Sens. 2022, 14, 488 3 of 23 method and the correlations with the influencing factors. Section 6 presents a discussion and summary of this work.

Instrumentation and Data
In this section, we introduce the theory of the spatial heterodyne spectroscopy (SHS) technique and the basic specifications of GMI, and we briefly outline the preprocessing of raw GMI data.

SHS Technology
SHS is a new interference spectroscopic technique developed by F. L. Roesler and J. Harlander [20]. When compared to a conventional Michelson interferometer, the SHS configuration replaces the return mirrors with diffraction gratings (G1 and G2, as shown in Figure 1). Incident light enters through aperture A1 and is collimated by lens L1; then, a wavelength-dependent Fizeau fringe pattern is produced by G1 and G2. Finally, this pattern is recorded on imaging detector I through the relay of lenses L2 and L3. 3 introduces the CBCR method. In Section 4, the influence of the surface bidirectional reflectance on the retrieval accuracy is analyzed from a theoretical perspective, the two methods are compared, and the main influencing factors are delineated. Section 5 compares the two retrieval methods using GMI measurements and analyzes the correction effect of the new method and the correlations with the influencing factors. Section 6 presents a discussion and summary of this work.

Instrumentation and Data
In this section, we introduce the theory of the spatial heterodyne spectroscopy (SHS) technique and the basic specifications of GMI, and we briefly outline the preprocessing of raw GMI data.

SHS Technology
SHS is a new interference spectroscopic technique developed by F. L. Roesler and J. Harlander [20]. When compared to a conventional Michelson interferometer, the SHS configuration replaces the return mirrors with diffraction gratings (G1 and G2, as shown in Figure 1). Incident light enters through aperture A1 and is collimated by lens L1; then, a wavelength-dependent Fizeau fringe pattern is produced by G1 and G2. Finally, this pattern is recorded on imaging detector I through the relay of lenses L2 and L3. The Fourier transformation of the fringe pattern recovers the spectrum. For the input spectrum ( ), the intensity on imaging detector I as a function of the position x is given by the following integral: where is the wavenumber of light, 0 is the Littrow wavenumber, and is the Littrow angle. A zero spatial frequency corresponds to the Littrow wavenumber of the gratings and does not correspond to a zero wavenumber, as is the approach in conventional Fourier transform spectroscopy (FTS).

GMI Specifications
GMI, which is based on SHS, boasts a high optical throughput and spectral resolution, and its structure is compact and stable. GMI observes solar light reflected from the The Fourier transformation of the fringe pattern recovers the spectrum. For the input spectrum B(σ), the intensity on imaging detector I as a function of the position x is given by the following integral: where σ is the wavenumber of light, σ 0 is the Littrow wavenumber, and θ is the Littrow angle. A zero spatial frequency corresponds to the Littrow wavenumber of the gratings and does not correspond to a zero wavenumber, as is the approach in conventional Fourier transform spectroscopy (FTS).

GMI Specifications
GMI, which is based on SHS, boasts a high optical throughput and spectral resolution, and its structure is compact and stable. GMI observes solar light reflected from the Earth's surface from near-infrared (NIR) wavelengths to the shortwave infrared (SWIR) region and measures raw interferograms, which are then converted into radiance spectra through a Fourier transformation. This instrument has one NIR band (centered at 0.76 µm, known as GMI band 1) and three SWIR bands (centred at 1.58, 1.65, and 2.0 µm, known as GMI bands Remote Sens. 2022, 14, 488 4 of 23 2, 3, and 4, respectively) with full width at half-maximum (FWHM) values of the instrument line shape (ILS) function of 0.035, 0.067, 0.067, and 0.113 nm, respectively ( Table 1).
The NIR band provides information on the surface pressure, clouds and aerosols. In addition, the SWIR-1 and SWIR-2 bands obtain important data regarding the weak absorption bands of CO 2 and CH 4 and provide information on CO 2 columns and CH 4 columns with high near-surface sensitivity. In contrast, the SWIR-3 band measures the strong absorption of CO 2 and accordingly provides additional CO 2 information.

Recovery of Spectral Interferograms
The raw interferograms measured by the SHS technique were converted into radiance spectra through quality filtering, error corrections, zero filling, a Fourier transformation, and spectral and absolute radiance calibrations. The detailed process is not described in this paper; rather, other publications by the team responsible for the GMI Level 1 product are responsible for thoroughly illustrating the preprocessing scheme [9]. After these preprocessing steps, the absolute radiance spectrum of each band was generated with 8192 spectral points, as shown in Figure 2. Earth's surface from near-infrared (NIR) wavelengths to the shortwave infrared (SWIR) region and measures raw interferograms, which are then converted into radiance spectra through a Fourier transformation. This instrument has one NIR band (centered at 0.76 μm, known as GMI band 1) and three SWIR bands (centred at 1.58, 1.65, and 2.0 μm, known as GMI bands 2, 3, and 4, respectively) with full width at half-maximum (FWHM) values of the instrument line shape (ILS) function of 0.035, 0.067, 0.067, and 0.113 nm, respectively (Table 1). The NIR band provides information on the surface pressure, clouds and aerosols. In addition, the SWIR-1 and SWIR-2 bands obtain important data regarding the weak absorption bands of CO2 and CH4 and provide information on CO2 columns and CH4 columns with high near-surface sensitivity. In contrast, the SWIR-3 band measures the strong absorption of CO2 and accordingly provides additional CO2 information.

Recovery of Spectral Interferograms
The raw interferograms measured by the SHS technique were converted into radiance spectra through quality filtering, error corrections, zero filling, a Fourier transformation, and spectral and absolute radiance calibrations. The detailed process is not described in this paper; rather, other publications by the team responsible for the GMI Level 1 product are responsible for thoroughly illustrating the preprocessing scheme [9]. After these preprocessing steps, the absolute radiance spectrum of each band was generated with 8192 spectral points, as shown in Figure 2.

Retrieval Method
The CBCR method was developed to retrieve the column-averaged dry-air mole fraction of carbon dioxide (XCO2) from the GMI radiance spectra in the NIR, SWIR-1, and SWIR-3 bands, which are characterized in Table 1. A flow chart of the CBCR method is

Retrieval Method
The CBCR method was developed to retrieve the column-averaged dry-air mole fraction of carbon dioxide (XCO 2 ) from the GMI radiance spectra in the NIR, SWIR-1, and SWIR-3 bands, which are characterized in Table 1. A flow chart of the CBCR method is depicted in Figure 3. The retrieval process begins with the inputs of high-quality GMI spectra filtered for retrievals. First, information about the GMI measurements and the corresponding a priori values of the state vectors from atmosphere and surface datasets are collected. Second, a forward model is used to generate synthetic radiance spectra and each corresponding Jacobian, which is the derivative of the radiance spectra, with respect to each state vector, for the GMI NIR, SWIR-1, and SWIR-3 bands with the recorded values. Finally, with the a priori values of the state vectors and the precalculated covariance matrix of the state vectors, a retrieval algorithm is applied to iteratively find the best solution of the state vectors. In each iteration, a forward model with new iterated inputs is invoked to generate new spectra until the convergence test is satisfied. Then, the retrieval process stops and outputs the retrieved XCO 2 . The remainder of this section describes each component of the retrieval algorithm.

Retrieval Method
The CBCR method was developed to retrieve the column-averaged dry-air mole fraction of carbon dioxide (XCO2) from the GMI radiance spectra in the NIR, SWIR-1, and SWIR-3 bands, which are characterized in Table 1. A flow chart of the CBCR method is depicted in Figure 3. The retrieval process begins with the inputs of high-quality GMI spectra filtered for retrievals. First, information about the GMI measurements and the corresponding a priori values of the state vectors from atmosphere and surface datasets are collected. Second, a forward model is used to generate synthetic radiance spectra and each corresponding Jacobian, which is the derivative of the radiance spectra, with respect to each state vector, for the GMI NIR, SWIR-1, and SWIR-3 bands with the recorded values. Finally, with the a priori values of the state vectors and the precalculated covariance matrix of the state vectors, a retrieval algorithm is applied to iteratively find the best solution of the state vectors. In each iteration, a forward model with new iterated inputs is invoked to generate new spectra until the convergence test is satisfied. Then, the retrieval process stops and outputs the retrieved XCO2. The remainder of this section describes each component of the retrieval algorithm.
GMI radiance spectra

Forward Model
The forward model was used to simulate GMI measurements and was constructed as shown in Figure 4. The main module of the GMI forward model is an atmospheric radiative transfer model (RTM), which simulates the effects of the atmosphere and ground surface on the solar radiance that propagates through the atmosphere and is reflected to

Forward Model
The forward model was used to simulate GMI measurements and was constructed as shown in Figure 4. The main module of the GMI forward model is an atmospheric radiative transfer model (RTM), which simulates the effects of the atmosphere and ground surface on the solar radiance that propagates through the atmosphere and is reflected to the instrument. Through the combination of a self-defined land surface model, a self-defined aerosol model, an instrument model and a high spectral resolution solar model in the RTM, the GMI forward model was constructed.  The RTM we adopted in this study was SCIATRAN v3.5, which was developed by V.V. Rozanov, and A.V. Rozanov from the Department of Remote Sensing, Institute of Environmental Physics, University of Bremen, Germany [21]. This model includes some common modules, such as the high-resolution transmission molecular absorption (HI-TRAN) database, solar spectrum database, atmospheric model, and surface model, and we could choose to input user-defined data. It also includes a line-by-line (LBL) gas absorption module. All of these modules are important for the simulation of GMI spectra. The RTM we adopted in this study was SCIATRAN v3.5, which was developed by V.V. Rozanov, and A.V. Rozanov from the Department of Remote Sensing, Institute of Environmental Physics, University of Bremen, Germany [21]. This model includes some common modules, such as the high-resolution transmission molecular absorption (HITRAN) database, solar spectrum database, atmospheric model, and surface model, and we could choose to input user-defined data. It also includes a line-by-line (LBL) gas absorption module. All of these modules are important for the simulation of GMI spectra.
In addition to the default solar models of SCIATRAN, we calculated a new highresolution solar spectrum for GMI through the multiplication of a solar continuum spectrum with an absorption spectrum. The former was calculated through the polynomial fitting of the low-resolution extraterrestrial solar spectrum supplied by Gérard Thuillier [22], while the latter (supplied by Geoff Toon) consisted of high-resolution solar pseudo-transmittance spectra in the range of 600-25,000 cm −1 [23].
We chose a kernel-driven bidirectional reflectance model to characterize the reflectance of the land surface. This kernel-driven bidirectional reflectance model can describe the BRDF of the land surface physically and is not complex, which is suitable for operational retrievals. The model given by Schaaf, shown as follows, is widely used [24]: where K iso is the isotropic scattering kernel; K vol is the volume scattering kernel; K geo is the geometric-optical scattering kernel; θ i , θ v , and ϕ are the solar zenith angle (SZA), viewing zenith angle (VZA), and relative azimuth angle, respectively; λ is the wavenumber; and f iso (λ), f vol (λ), and f geo (λ) are the kernel coefficients, which represent the weights of isotropic scattering, volume scattering, and geometric-optical scattering, respectively. In the above formula, the isotropic scattering, volume scattering and geometric-optical scattering weights are the coefficients of the corresponding three kernels and are functions related to the band. The maximum value of any weight is 1. The common RossThick (Ross) kernel and LiSparseR (Li) kernel functions are combined into the Ross-Li model, where the Ross kernel function is expressed as: where g is the scattering angle, cos g = cos θ i cos θ v + sin θ i sin θ v cos ϕ. The Li kernel functions are expressed as Therein: where the ratio of h to b is the height-to-width ratio of the canopy ellipsoid, which is generally fixed at 2.0. We defined 11 aerosol types in combination with the 550 nm aerosol optical depth (AOD) to characterize aerosols. The global aerosol distribution is relatively complex in time and space. Considering the feasibility of acquiring the appropriate data sources, the statistics of the aerosol types and their AOD distribution characteristics were obtained through global aerosol cluster analysis. Taylor et al. [25] obtained a global aerosol mixing model through the cluster analysis of 7 years of Goddard Chemistry Aerosol Radiation and Transport (GOCART) model data. The GOCART model [26] divides global aerosols into sulfate (Su) aerosols, black carbon (BC) aerosols, organic carbon (OC) aerosols, dust (DU) aerosols, and sea salt (SS), and BC aerosols and OC aerosols can be combined into biomassburning (BB) aerosols. We considered the temporal variation of the aerosol distribution according to the four seasons: March, April, and May were spring; June, July, and August were summer; September, October, and November were autumn; and December, January, and February were winter. The data were reclustered to obtain the distributions of the aerosol types in all four seasons around the world, as shown in Figure 5. The compositions of the various seasonal aerosol models are shown in Table 2. Prior knowledge of the global distributions of these aerosol models provides support for the global GMI retrieval of CO 2 .
where the ratio of ℎ to is the height-to-width ratio of the canopy ellipsoid, which is generally fixed at 2.0.
We defined 11 aerosol types in combination with the 550 nm aerosol optical depth (AOD) to characterize aerosols. The global aerosol distribution is relatively complex in time and space. Considering the feasibility of acquiring the appropriate data sources, the statistics of the aerosol types and their AOD distribution characteristics were obtained through global aerosol cluster analysis. Taylor et al. [25] obtained a global aerosol mixing model through the cluster analysis of 7 years of Goddard Chemistry Aerosol Radiation and Transport (GOCART) model data. The GOCART model [26] divides global aerosols into sulfate (Su) aerosols, black carbon (BC) aerosols, organic carbon (OC) aerosols, dust (DU) aerosols, and sea salt (SS), and BC aerosols and OC aerosols can be combined into biomass-burning (BB) aerosols. We considered the temporal variation of the aerosol distribution according to the four seasons: March, April, and May were spring; June, July, and August were summer; September, October, and November were autumn; and December, January, and February were winter. The data were reclustered to obtain the distributions of the aerosol types in all four seasons around the world, as shown in Figure 5. The compositions of the various seasonal aerosol models are shown in Table 2. Prior knowledge of the global distributions of these aerosol models provides support for the global GMI retrieval of CO2.   The instrument model consisted of an ILS model and a noise model. The simulated hyper-resolution radiance spectrum was convolved with the ILS of the corresponding GMI band, whereas the instrument noise model contained instrument noise and forward model errors. In this study, we neglected forward model errors and took only instrument noise into account for simplicity. Before the launch of GMI, a careful set of calibration tests were performed to measure the ILS and instrument noise, which were used directly for this study, under the assumption that they remained unchanged in the first half-year of on-orbit operation [9].

State Vectors and Auxiliary Parameters
State vectors and auxiliary parameters all influence the radiance spectra; the former represents the parameters that are optimized and retrieved in the inversion, while the auxiliary parameters are not. The state vectors and auxiliary parameters employed in the proposed GMI CBCR method are shown in Table 2. The state vectors included the CO 2 profile, three kernel coefficients of the SWIR-1 band, albedos of the NIR and SWIR-3 bands, 550 nm AOD and wavenumber factors. Only the BRDF model was adopted for the SWIR-1 band because of the extremely time-consuming radiation calculation, and the manually calculated coefficients of the three kernel functions are functions of synchronous calculations that SCIATRAN does not currently have. The auxiliary parameters included the aerosol type, surface height, temperature profile, H 2 O profile, and pressure profile.
The a priori CO 2 profiles were derived from the Goddard Earth Observing System (GEOS)-Chem model [27], which is a global 3D chemical transport model driven by GEOSassimilated meteorological fields from the NASA Global Modelling and Assimilation Office [28]. GEOS-Chem simulates CO 2 at a 2 • × 2.5 • (latitude × longitude) horizontal resolution and 47 vertical pressure levels from 1005.65 hPa to 0.038 hPa. We divided the atmosphere into 46 levels from the ground surface to the top of the atmosphere, according to the standard US76 atmospheric model, and interpolated the simulated CO 2 profiles into these layers. To acquire the a priori 2020 profiles, considering that the supportive datasets for GEOS-Chem are delayed by approximately two years, we used GEOS-Chem to calculate CO 2 over five years from 2015 to 2019 and then extrapolated the values to 2020.
The a priori covariance matrix of CO 2 profiles was defined as a function of the seasonal bias plus the interannual variance and the synoptic covariance matrix [29], and each element was calculated at each of the 46 vertical levels using five years of GEOS-Chem CO 2 results. The covariance varied greatly with altitude; the total a priori uncertainty was the largest near the surface at 94 ppm 2 and decreased with increasing altitude, reaching a minimum of 1 ppm 2 . We multiplied the calculated a priori uncertainties of the CO 2 profiles by a factor of 10 to avoid excessively strong constraints.
The parameters of the three kernel functions in the SWIR-1 band were derived from the MODIS 16 day synthetic grid product MCD43, which provided the Ross-Li model coefficients in MODIS bands 1-7. According to the spatial resolution of the product, the data were divided into the MCD43A1 product with a 0.5 km resolution and the MCD43C1 product with a 0.05 • resolution. This paper used band 6 (1.628-1.652 µm) of the MCD43A1 product with a relatively fine resolution for the GMI SWIR-1 band. We selected a typical urban area and a typical vegetation area in Beijing, China, as an example in Figure 6a to analyze the BRDF coefficient characteristics of different underlying surface types. Point A was in the central area of Beijing, and point B was in a suburban vegetation area for comparison. The GMI observation points are projected as red circles with a diameter of 10 km on the ground, so it was necessary to resample the MCD43A1 data with a resolution of 0.5 km. The blank areas in the red circles represent the lack of data within the product, mainly due to cloud cover and snow cover; hence, the AOD was too large and could not be retrieved. In addition, because snow falls in Beijing during winter, snow pixel preprocessing was needed. due to cloud cover and snow cover; hence, the AOD was too large and could not be retrieved. In addition, because snow falls in Beijing during winter, snow pixel preprocessing was needed. Figure 6b shows the variations in the three core BRDF coefficients from 2012 to 2016. The left one corresponds to the urban area and the right one corresponds to the vegetation area. Although the three core BRDF coefficients in the urban area changed seasonally, the variation range is small; in contrast, in the vegetation area, the coefficient values vary greatly with the season. When compared with that in vegetation areas, the underlying surface structure in urban areas does not change greatly with the seasons, resulting in small variations in the volume scattering and geometric-optical scattering coefficients, and thus, their change trends are the same as that of the isotropic scattering coefficient. As a reference, the vegetation area does not exhibit the abovementioned variations due to the seasonal changes in vegetation cover. Therefore, priori values in winter and spring are different from summer and autumn for vegetation regions.  The covariance of each BRDF kernel function coefficient was calculated by performing a statistical analysis of the temporal variation of the coefficient distribution in the last five years. For regions in which the BRDF coefficients exhibited seasonal characteristics, the a priori covariance was calculated according to winter-spring and summer-autumn.  Figure 6b shows the variations in the three core BRDF coefficients from 2012 to 2016. The left one corresponds to the urban area and the right one corresponds to the vegetation area. Although the three core BRDF coefficients in the urban area changed seasonally, the variation range is small; in contrast, in the vegetation area, the coefficient values vary greatly with the season. When compared with that in vegetation areas, the underlying surface structure in urban areas does not change greatly with the seasons, resulting in small variations in the volume scattering and geometric-optical scattering coefficients, and thus, their change trends are the same as that of the isotropic scattering coefficient. As a reference, the vegetation area does not exhibit the abovementioned variations due to the seasonal changes in vegetation cover. Therefore, priori values in winter and spring are different from summer and autumn for vegetation regions.
The covariance of each BRDF kernel function coefficient was calculated by performing a statistical analysis of the temporal variation of the coefficient distribution in the last five years. For regions in which the BRDF coefficients exhibited seasonal characteristics, the a priori covariance was calculated according to winter-spring and summer-autumn. For satellite observations with a spatial resolution of 10 km, the change in the surface structure is often small, and therefore, the corresponding BRDF kernel function coefficient is also small. As shown in Figure 6, the standard deviations of f iso , f vol , and f geo in the urban area are small, at 0.02, 0.018 and 0.012, respectively; whereas, the corresponding standard deviations in the vegetation area are large, at 0.032, 0.021 and 0.013 in winter-spring, and 0.02 in summer-autumn.
The initial reflectances in the NIR and SWIR-3 bands originated from the MODIS MCD43C3 product, which has a spatial resolution of 0.05 • and a temporal resolution of 16 days for all 7 bands. In this paper, the Band 2 (NIR, 0.841-0.876 µm) and Band 7 (SWIR-3, 2.105-2.155 µm) products were used. Using the statistics of the temporal variations of the MCD43C3 product data over the last three years, the reflectance values corresponding to the observation times were obtained by fitting, and the element values of the a priori covariance were calculated according to the covariance matrix.
The prior AOD values were derived from the monthly average aerosol product of MODIS MOD08, which has a spatial resolution of 1 km and can provide the 550 nm AOD. The 550 nm AOD values corresponding to the observation times were obtained by fitting the time-varying statistics of the data over the last three years, and the a priori covariance was calculated according to the covariance matrix. The variation in AOD with time in some areas is often complicated. Therefore, the averaged fitting of MODIS monthly products with different orbits was obtained. Considering that DPC aerosol products are not routinely produced at present, it will be helpful to improve the retrieval ability by adopting the current method and adopting a higher accuracy of AOD data sources in the future.
Wavenumber factors were also included among the state vectors because spectral shifts are unavoidable in satellite spectral measurements. The a priori wavenumber shift and squeeze factor were established according to the method developed by Ye [30], and the a priori uncertainties were set to 1. The second-order polynomial coefficients were retrieved simultaneously. The a priori wavenumber shift and squeeze factor were both set to 0, and the a priori uncertainties were set to 1.
Aerosol types were auxiliary data, and the initial values were obtained from a global aerosol type lookup table according to the seasons and spatial locations of the observation points. The spatial resolution of this lookup table, which was divided into four seasons, was 1 km. See Section 3.1 for further details.
To acquire the elevations at the observation points, we adopted Shuttle Radar Topography Mission (SRTM) surface elevation data with a spatial resolution of 3 arc seconds (approximately 90 m) [31]. According to the spatial positions of the observation points, all the elevation data in the field of view within a diameter of 10 km centered on the observation point were obtained, and the corresponding elevations of the observation points were obtained through averaging.
The temperature profile, H 2 O profile, and pressure profile are from the European Centre for Medium-Range Weather Forecasts (ECMWF) (before March 2019) and ERA5 (after March 2019). The spatial resolution of these reanalysis data was 0.125 • , and the vertical stratification comprised 37 layers ranging from 1 to 1000 hPa. The upper atmospheric profile of the atmospheric model, corresponding to the integration of latitude zone and season, is supplemented. In this way, the completed profile for each observation point from the surface to the top of the atmosphere was obtained.

Retrieval Algorithm
GMI measurement y can be expressed as follows: where y contains the radiance spectra of the three GMI bands (NIR, SWIR-1, and SWIR-3), the forward model F is used to simulate the GMI measurements, illustrated in Section 3.1, x is the state vector, b represents the auxiliary parameters, with the elements illustrated in Section 3.2, and ε contains the instrument noise and forward model errors. The instrument noise for GMI was derived from laboratory measurements performed before launching the instrument; hence, this noise was considered to remain unchanged during the first half-year of on-orbit operation [9]. The optimal estimation method retrieves the solution of state vectors (x) that are optimized to yield simulated spectra that best match the measured spectra, and are simultaneously constrained by a priori information with a suitable error covariance [32]. First, the cost function χ 2 , which represents the difference between the measured spectra y and the simulated spectra F(x), as well as the difference between the retrieved state vector x and the a priori state vector x a , is constructed as follows: where S ε is the measurement error covariance matrix, S a is the a priori covariance matrix, and the superscripts T and −1 represent matrix transpose and matrix retrieval operations, respectively. The solution of the state vector x minimizes the cost function with the maximum a posteriori probability. The Levenberg-Marquardt method is used to update the state vector using the following iterative equation at each step: where the subscript i represents the i-th iteration, K = ∂F(x)/∂x is the Jacobian representing the first derivative of F(x), and γ is the Levenberg-Marquardt parameter. The value of γ is initially set to 10.0 and is updated according to the strategy proposed by C.D. Rodgers based on the ratio R of the actual reduction in χ 2 to the forecasted reduction in χ 2 under the assumption of linearity; this ratio is expressed as follows: During the iteration, if the updated value of the state vector was less than 10 −4 of the a priori value, the iteration was considered convergent, and we obtained the retrievals successfully.
The retrieved XCO 2 profilesx are multiplied by the pressure profile P assigned to the atmospheric level to calculate the dry-air mixing ratio X CO 2 using the following equation:

Simulations of GMI Observations
To quantitatively study the influence of atmospheric scattering corresponding to the surface bidirectional reflectance characteristics on the CO 2 retrieval accuracy, GMI data were first simulated. The urban area in Figure 6a is taken as an example. Urban areas often experience various changes in their atmospheric conditions, and the simulation (as shown in Figure 7) shows that this urban area generally exhibits remarkable bidirectional reflectance characteristics. Assuming that the initial concentration of CO 2 is 400 ppm, the surface bidirectional reflectance characteristics adopt the BRDF parameter values of urban areas provided by MODIS, as shown in Figure 6b

Retrieval Results
The GMI_XCO2 method regards the surface reflectance and the corresponding atmospheric scattering as slow components and reduces the influence of atmospheric scattering by removing these slow components from the logarithmic radiance spectrum. In contrast, the CBCR method adds the surface BRDF model to the forward model in the retrieval process and simultaneously inverts the core BRDF coefficients and AOD of the surface BRDF model in the retrieval of the CO2 content to correct for both the surface bidirectional reflectance characteristics and the influence of the corresponding atmospheric scattering. Under an AOD of 0.05, the correction results of the GMI_XCO2 method are shown as the curves with black squares in Figure 8. The influence of the surface bidirectional reflectance on the CO2 retrieval accuracy has certain directional characteristics, but the error is small, basically within 0.6 ppm, and the maximum error is less than 1 ppm. As shown by the curves with red circles in Figure 8, the retrieval results obtained by the CBCR method are very stable under various observation conditions; the retrieval errors caused by changes in the surface reflectivity with the SZA and VZA are well corrected, and the errors are all less than 0.3 ppm.  Figure 7 shows the corresponding BRDF characteristics of the underlying urban surface of central Beijing in all four seasons with varying SZA. The reflectance characteristics of the surface obviously change with the SZA and VZA, especially the difference in reflectivity corresponding to the forward observation on the principal plane, which is the most significant. To ascertain the influence of the surface bidirectional reflectance characteristics on atmospheric CO 2 retrievals, quantitative research was carried out under different SZAs and VZAs based on the principal plane.

Retrieval Results
The GMI_XCO 2 method regards the surface reflectance and the corresponding atmospheric scattering as slow components and reduces the influence of atmospheric scattering by removing these slow components from the logarithmic radiance spectrum. In contrast, the CBCR method adds the surface BRDF model to the forward model in the retrieval process and simultaneously inverts the core BRDF coefficients and AOD of the surface BRDF model in the retrieval of the CO 2 content to correct for both the surface bidirectional reflectance characteristics and the influence of the corresponding atmospheric scattering. Under an AOD of 0.05, the correction results of the GMI_XCO 2 method are shown as the curves with black squares in Figure 8. The influence of the surface bidirectional reflectance on the CO 2 retrieval accuracy has certain directional characteristics, but the error is small, basically within 0.6 ppm, and the maximum error is less than 1 ppm. As shown by the curves with red circles in Figure 8, the retrieval results obtained by the CBCR method are very stable under various observation conditions; the retrieval errors caused by changes in the surface reflectivity with the SZA and VZA are well corrected, and the errors are all less than 0.3 ppm. When aerosol scattering was neglected or the effect of aerosol scattering was small, the retrieval error caused by the surface bidirectional characteristics is relatively minor. However, if only the AOD was considered in the CO2 retrieval process, the data utilization efficiency is greatly reduced. Accordingly, we assumed that when the AOD is 0.4, other conditions are set in the same way as before, and the statistics of the XCO2 retrieval errors corresponding to the GMI_XCO2 and CBCR methods are shown in Figure 9. As shown in Figure 9a, in the backward observations, the CO2 retrieval error caused by atmospheric scattering, corresponding to the two-way surface reflectance characteristics, increases with increasing SZA. With increasing SZA, the length of the optical absorption path of CO2 increases, which is more affected by aerosol scattering, with an average value of −4.91 ppm. In the forward observations, with increasing SZA, the retrieval error first decreases and then increases. With increasing VZA, the retrieval error first decreases and then increases rapidly, and the maximum retrieval error reaches 20 ppm. Figure 9b shows the retrieval errors after correcting for the surface bidirectional reflectance characteristics and atmospheric scattering. When the VZA is less than 50°, the retrieval result is stable, and When aerosol scattering was neglected or the effect of aerosol scattering was small, the retrieval error caused by the surface bidirectional characteristics is relatively minor. However, if only the AOD was considered in the CO 2 retrieval process, the data utilization efficiency is greatly reduced. Accordingly, we assumed that when the AOD is 0.4, other conditions are set in the same way as before, and the statistics of the XCO 2 retrieval errors corresponding to the GMI_XCO 2 and CBCR methods are shown in Figure 9. As shown in Figure 9a, in the backward observations, the CO 2 retrieval error caused by atmospheric scattering, corresponding to the two-way surface reflectance characteristics, increases with increasing SZA. With increasing SZA, the length of the optical absorption path of CO 2 increases, which is more affected by aerosol scattering, with an average value of −4.91 ppm.
In the forward observations, with increasing SZA, the retrieval error first decreases and then increases. With increasing VZA, the retrieval error first decreases and then increases rapidly, and the maximum retrieval error reaches 20 ppm. Figure 9b shows the retrieval errors after correcting for the surface bidirectional reflectance characteristics and atmospheric scattering. When the VZA is less than 50 • , the retrieval result is stable, and the error hardly changes with a change in the SZA. The correction effect is optimal when the SZA is less than 30 • , and the retrieval error is stable, with an average retrieval error of −1.4 ppm; when the VZA is greater than 60 • , the error increases to a maximum of 2.2 ppm. When the SZA is greater than 30 • , the retrieval result is relatively stable, but the retrieval error increases. When the SZA is 60 • , the average retrieval error is 2.98 ppm, and the maximum retrieval error is 5.2 ppm when the VZA is 70 • .

Data Preparation
GMI started to routinely collect global observations in August 2018. In the early stage, we processed the Level One data and obtained stable spectral data. As a result, the spectral data from August 2018 to April 2020 were completely processed and verified. Two screening steps were implemented before the GMI L1 product inversion. The first step was to screen the time, location, and SZA; land observations with SZA < 70° were selected from August 2018 to April 2020. The second step was to eliminate data contaminated by cloud cover. At present, the GMI NIR spectrum is used to monitor clouds combined with DPC cloud monitoring on the same platform [33]. Considering that both monitoring methods are immature and that possible cloud cover can interfere with the evaluation and analysis of the CO2 product quality, at present, any retrieval judged as being influenced by clouds was eliminated.
The ground-based FTS XCO2 measurements from Total Carbon Column Observing Network (TCCON) were used here to validate the GMI retrievals because these reference data show the highest accuracy of ~0.8 ppm [34] and therefore provide a uniform method to analyze all the stations [35]. Therefore, TCCON data (https://doi.org/10.14291/TCCON.GGG2014, last accessed 1 November 2021) were used to evaluate the GMI retrieval results. The GMI XCO2 retrievals within the square lati-

Data Preparation
GMI started to routinely collect global observations in August 2018. In the early stage, we processed the Level One data and obtained stable spectral data. As a result, the spectral data from August 2018 to April 2020 were completely processed and verified. Two screening steps were implemented before the GMI L1 product inversion. The first step was to screen the time, location, and SZA; land observations with SZA < 70 • were selected from August 2018 to April 2020. The second step was to eliminate data contaminated by cloud cover. At present, the GMI NIR spectrum is used to monitor clouds combined with DPC cloud monitoring on the same platform [33]. Considering that both monitoring methods are immature and that possible cloud cover can interfere with the evaluation and analysis of the CO 2 product quality, at present, any retrieval judged as being influenced by clouds was eliminated.
The ground-based FTS XCO 2 measurements from Total Carbon Column Observing Network (TCCON) were used here to validate the GMI retrievals because these reference data show the highest accuracy of~0.8 ppm [34] and therefore provide a uniform method to analyze all the stations [35]. Therefore, TCCON data (https://doi.org/10.14291/TCCON. GGG2014, last accessed 1 November 2021) were used to evaluate the GMI retrieval results. The GMI XCO 2 retrievals within the square latitude/longitude regions centered on the TCCON stations and spaced at ±5 • were adopted, and the TCCON data within 1 h of the GMI overpass time (approximately 13:30 local time) were averaged as a reference. From August 2018 to April 2020, 12 TCCON stations with a large amount of data were available, as shown in Table 3.

Retrieval Results and Validation
For the GMI observation data from the 12 TCCON sites, the XCO 2 results obtained by the CBCR method and the GMI_XCO 2 method are compared in Figure 10. The yellow points are the XCO 2 results of the CBCR method, the green points are the GMI_XCO 2 retrieval results, the blue points are all the TCCON observations from August 2018 to April 2020, and the red line in each panel represents the temporal variation trend of atmospheric CO 2 obtained by fitting all of the TCCON stations' observations from the beginning date.

Retrieval Results and Validation
For the GMI observation data from the 12 TCCON sites, the XCO2 results obtained by the CBCR method and the GMI_XCO2 method are compared in Figure 10. The yellow points are the XCO2 results of the CBCR method, the green points are the GMI_XCO2 retrieval results, the blue points are all the TCCON observations from August 2018 to April 2020, and the red line in each panel represents the temporal variation trend of atmospheric CO2 obtained by fitting all of the TCCON stations' observations from the beginning date.  Figure 10 shows that the atmospheric CO2 retrieval results obtained by the CBCR method at all stations exhibited significant seasonal variation characteristics and annual growth trends that are consistent with the temporal variation characteristics of atmospheric CO2 reflected by the TCCON station measurements; the lowest correlation is 45.75%, and the highest correlation is 76.8%, as shown in Table 4. Since there were few points of consistent time between CBCR method retrievals and TCCON measurements, especially in winter, there are large vacancies in CBCR, so the time fitting lines of the two are quite different, which caused low correlation coefficients. In contrast, the atmospheric CO2 retrieval results obtained by the GMI_XCO2 method are discrete, yielding low correlations with the temporal variation characteristics of atmospheric CO2 ranging from 17.05 to 37.4%; these low correlations make it difficult to reveal the clear seasonal variations and annual growth characteristics of atmospheric CO2. The CBCR and GMI_XCO2 method inversions were accurately verified using TCCON ground measurements corresponding to the GMI satellite observations. The   Figure 10 shows that the atmospheric CO 2 retrieval results obtained by the CBCR method at all stations exhibited significant seasonal variation characteristics and annual growth trends that are consistent with the temporal variation characteristics of atmospheric CO 2 reflected by the TCCON station measurements; the lowest correlation is 45.75%, and the highest correlation is 76.8%, as shown in Table 4. Since there were few points of consistent time between CBCR method retrievals and TCCON measurements, especially in winter, there are large vacancies in CBCR, so the time fitting lines of the two are quite different, which caused low correlation coefficients. In contrast, the atmospheric CO 2 retrieval results obtained by the GMI_XCO 2 method are discrete, yielding low correlations with the temporal variation characteristics of atmospheric CO 2 ranging from 17.05 to 37.4%; these low correlations make it difficult to reveal the clear seasonal variations and annual growth characteristics of atmospheric CO 2 .
The CBCR and GMI_XCO 2 method inversions were accurately verified using TCCON ground measurements corresponding to the GMI satellite observations. The mean deviations and standard deviations between the two algorithms' retrieval results and the TCCON measurements are shown in Table 5. The calculated values reveal bias in the atmospheric CO 2 retrievals of the GMI_XCO 2 method: all sites had a mean deviation of 0.58 ppm, and the standard deviation was 5.64 ppm, which far exceeded the accuracy requirements of satellite retrievals (4 ppm). Furthermore, different sites had a large difference (minimum of 3.58 ppm and maximum of 6.2 ppm, with a difference of 2.62 ppm). In addition, the retrieval results of the CBCR method also had an obvious negative bias: the mean deviation among all sites was −1.33 ppm, but the standard deviation was reduced to 3.13 ppm. Moreover, the retrieval error was consistent among the stations with the CBCR method; the largest retrieval error was 3.30 ppm, the lowest retrieval error was 2.43 ppm, and the maximum difference was only 0.87 ppm. Likewise, the standard deviations in the atmospheric CO 2 retrieval results of the CBCR and GMI_XCO 2 methods under different AODs, surface reflectances, SZAs, and VZAs were also compared, as shown in Figure 11. According to the distribution ranges of the AODs, surface reflectivity, SZAs, and VZAs corresponding to the GMI observation points, the AOD was set from 0 to 0.5, with an interval of 0.1 (very few AODs exceed 0.5, and all of the remainder are below 0.5). The surface reflectivity was similarly set from 0 to 0.5 with an interval of 0.1. The SZA ranged from 20 • to 70 • with an interval of 10 • . According to the GMI characteristics, the VZA was set to one of three values, namely, 3, 18, or 36 • , and adjacent values were merged into one category. The standard deviation was calculated according to a comparison between the GMI survey point retrieval results and the actual TCCON measurements within the interval corresponding to each element, and the mean deviation of each station was subtracted. also tends to increase with increasing AOD (with a slope of 3.082), but when compared with the GMI_XCO2 method, the proposed method greatly reduces the correlation between the deviation of the retrieval results and the AOD and significantly reduces the influence of the SZA on the retrieval results, with a slope of 0.002. Ultimately, the CBCR method significantly reduces the influences of the surface bidirectional reflectance characteristics under high AOD and high SZA conditions and improves the retrieval accuracy at each station. Figure 11. Correlation of the atmospheric CO2 retrieval deviation with AOD (a), SZA (b), viewing zenith angle (VZA) (c), and reflectance (d).

Conclusions
To retrieve high-precision CO2 products from GMI measurements and correct for the effect of atmospheric scattering induced by the surface BRDF, we developed a CO2 retrieval method, namely, the CBCR method, which retrieves CO2 together with the BRDF coefficients in the SWIR-1 band, the albedo in the NIR and SWIR-3 bands, the 550 nm AOD, and wavenumber factors simultaneously. The proposed method was coupled with Figure 11. Correlation of the atmospheric CO 2 retrieval deviation with AOD (a), SZA (b), viewing zenith angle (VZA) (c), and reflectance (d).
The standard deviations (STD) of the GMI_XCO 2 method retrieval results are shown as the blue points in Figure 11, while those of the CBCR method retrieval results are shown as the red triangles. The deviations of the GMI_XCO 2 method retrieval results display a certain correlation with the AOD and SZA. As the main trend, the larger the AOD is, the stronger the deviation of the inversion, and the slope with the AOD is 5.846. Similarly, the larger the SZA, the stronger the deviation of the inversion becomes, but the slope with the SZA is only 0.091. At the same time, the correlations between the SZA and the single surface reflectivity and VZA are not strong. The retrieval errors corresponding to different surface reflectances and different VZAs are basically the same, which is mainly due to the degree of influence of the surface bidirectional reflectance characteristics, which is closely related to the atmospheric conditions and is consistent with the simulation results under low and strong aerosol conditions. The deviation of the CBCR method retrieval results also tends to increase with increasing AOD (with a slope of 3.082), but when compared with the GMI_XCO 2 method, the proposed method greatly reduces the correlation between the deviation of the retrieval results and the AOD and significantly reduces the influence of the SZA on the retrieval results, with a slope of 0.002. Ultimately, the CBCR method significantly reduces the influences of the surface bidirectional reflectance characteristics under high AOD and high SZA conditions and improves the retrieval accuracy at each station.

Conclusions
To retrieve high-precision CO 2 products from GMI measurements and correct for the effect of atmospheric scattering induced by the surface BRDF, we developed a CO 2 retrieval method, namely, the CBCR method, which retrieves CO 2 together with the BRDF coefficients in the SWIR-1 band, the albedo in the NIR and SWIR-3 bands, the 550 nm AOD, and wavenumber factors simultaneously. The proposed method was coupled with the Ross-Li bidirectional reflectance model to characterize the reflectance of the land surface, and 11 seasonal aerosol types were combined with the 550 nm AOD to characterize aerosols. At the same time, the statistical analysis of satellite data was performed to acquire rich a priori information and to achieve the stable inversion of multiple parameters.
We examined the effect of atmospheric scattering due to the surface BRDF on the atmospheric CO 2 retrievals. Urban areas with obvious BRDF features were selected, and different AOD conditions (as low as 0.05 and as high as 0.4) were set. Several typical geometric conditions were investigated. The retrieval errors of the GMI_XCO 2 method (without considering the influence of atmospheric scattering corresponding to surface bidirectional reflectance characteristics) and the CBCR method were compared on the principal plane with obvious BRDF characteristics. The outcome of this comparison shows that at an AOD of 0.05, the influence of the BRDF characteristics on the CO 2 inversion exhibits a certain directionality, basically within 0.6 ppm; the maximum error is less than 1 ppm, and the error of the CBCR method after improvement is controlled within 0.3 ppm. However, under an AOD of 0.4, the atmospheric scattering influence corresponding to the surface bidirectional reflectance characteristics increased with increasing SZA in the backward observations, and the average CO 2 retrieval error was −4.91 ppm. In the forward observations, with increases in the SZA and VZA, the atmospheric scattering influence first decreased and then increased, and the maximum retrieval error could reach 20 ppm. Nevertheless, under this condition, when the SZA was less than 60 • and the VZA was less than 50 • , the retrieval results of the CBCR method were consistent, and the maximum retrieval error was less than 2.2 ppm, indicating that the proposed method can significantly improve the aerosol and geometric correlations for the CO 2 retrieval error caused by atmospheric scattering attributable to the surface bidirectional reflectance characteristics.
We used the GMI_XCO 2 and CBCR methods to invert GMI-measured spectral data from August 2018 to April 2020 and analyzed the retrieval results by using a large amount of measured data from 12 TCCON stations during the same period. The correlation of the atmospheric CO 2 retrieved by the GMI_XCO 2 method with the temporal variation of actual atmospheric CO 2 was as low as 34.7%, while the correlation increased to 76.8% for the CBCR method. There is a certain positive systematic bias in the atmospheric CO 2 of the GMI_XCO 2 method. The mean deviation of all stations was 0.58 ppm, and the standard deviation was as high as 5.64 ppm, which far exceeded the accuracy requirement (4 ppm) of satellite inversion. The mean deviation of different stations was quite different, ranging from 3.58 ppm to 6.2 ppm, and the difference was as high as 2.62 ppm. In contrast, there was a negative systematic bias in the retrieval results of the CBCR method, with a mean deviation of −1.33 ppm and an average standard deviation of 3.13 ppm at all stations. The retrieval error was consistent among the stations, with a maximum retrieval error of 3.3 ppm, a minimum retrieval error of 2.43 ppm, and a maximum difference of only 0.87 ppm.
By analyzing the variations in the deviations of the retrieval results of the GMI_XCO 2 and CBCR method with the AOD, surface reflectivity, SZA and VZA, we found that the retrieval error of the GMI_XCO 2 method exhibits certain correlations with the AOD and SZA: with increases in the AOD and SZA, the retrieval deviation became stronger, with the slopes corresponding to the AOD and SZA being 5.846 and 0.091, respectively. At the same time, the correlations with the single-surface reflectivity and VZA were not strong, mainly because the influence of the surface bidirectional reflectance characteristics on the atmospheric CO 2 retrievals is closely related to the aerosol and geometric conditions. The CBCR method reduced the correlations between the deviation of the CO 2 retrieval results and the AOD and SZA, with the slopes corresponding to the AOD and SZA being 3.082 and 0.002, respectively, suggesting that the proposed method significantly reduces the influence of surface bidirectional reflectance characteristics under high AOD and high SZA conditions. Ultimately, the CBCR method improves the ability to correct the atmospheric scattering influence due to surface bidirectional reflectance characteristics, provides an effective method to establish a high-precision global product for GMI measurements and applications, and demonstrates that describing the surface reflectance characteristics by the BRDF is a promising idea in the field of satellite CO 2 retrievals. In future work, we will improve the problems that remain in the CBCR method, such as its extremely time-consuming calculations and the strong dependence of the inversion on a priori constraints, to expand the inversion of global observations and to obtain more satellite observation data that are beneficial to the quantification of global CO 2 sources and sinks.