2.1. PICSs Overview
The Committee on Earth Observation Satellites (CEOS) recommends six PICSs for the long-term on-orbit radiometric calibration of Earth observation optical sensors. These sites possess several favorable characteristics, including temporal stability, spatial uniformity, and homogeneity [
13,
14].
One of these PICSs is Libya 4, which is located in the Libyan Desert, with consistently high reflectance throughout the year. Libya 4 is widely recognized as one of the most stable PICSs and has been extensively utilized [
11,
15]. Following CEOS recommendations, this study considered a circular area with a center at coordinates (28.55°N, 23.39°E) and a diameter of 20 km [
13].
Figure 1 is an image generated using the near-infrared (NIR) band (865 nm) of Landsat 8’s TOA reflectance that was captured on 20 August 2018. The region of interest is highlighted in red.
Figure 2 shows the stability of the TOA reflectance in the Libyan Desert from 1985 to 2021. The TOA reflectance data in the red band are extracted from the Landsat series of satellites. There are minor variations in the relative spectral response profiles among corresponding bands from different Landsat sensors [
16]. The data from Landsats 4, 5, and 7 were spectrally adjusted to match Landsat 8 OLI band 4 [
17]. The TOA reflectance does not exhibit any significant trend.
2.2. Data Sources
The study presented in this article relied on data from the following five satellite sensors: Terra/Aqua MODIS, Landsat 8 OLI, Suomi NPP VIIRS, and EO-1 Hyperion. The spatial resolutions and temporal ranges of the sensors used in this study are presented in
Table 1.
MODIS is a remote sensing sensor developed by NASA for studying the Earth’s atmospheric system. This study used MODIS data from the Terra (morning overpass) and Aqua (afternoon overpass) satellites, which image the entire Earth every one to two days [
18,
19].
MODIS has a radiometric resolution of 12 bits and offers up to 36 spectral bands [
20,
21]. The data used in this study for establishing the TOA reflectance model are derived from the first seven bands, covering the visible to SWIR region. The spatial resolution of these bands is either 250 m (bands 1 and 2) or 500 m (bands 3–7), as shown in
Table 2.
The MODIS data products used in this study include version 6.1 of MOD021KM, MYD021KM, and MYD08_D3. These products resample the data from bands 1 to 7 to a spatial resolution of 1000 m. MOD021KM/MYD021KM provides level 1B TOA radiance, reflectance, and related parameters [
22,
23]. MYD08_D3 provides daily atmospheric data, including AOD, water vapor, and ozone, which are used for atmospheric correction of the TOA reflectance model [
24].
Figure 3 illustrates the long-term TOA reflectance for bands 1 to 7 of Aqua MODIS over the Libyan Desert from 2003 to 2012. The TOA reflectance remained stable over time for each band.
Landsat 8 was launched in 2013 and is currently operated by the United States Geological Survey (USGS). It carries the OLI instrument, which has a swath width of 185 km and a length of 180 km, with a ground resolution of 30 m [
25]. This study focused on the first seven spectral bands for validation, covering wavelengths from visible to SWIR, as shown in
Table 3.
Suomi NPP is a low-Earth-orbit meteorological satellite that was launched in 2011. One of the key instruments on board Suomi NPP is VIIRS, which is a cross-track scanning radiometer. It provides twice-daily global coverage of the Earth’s surface and collects imagery data at a moderate spatial resolution of 750 m. VIIRS has 22 spectral bands, covering the wavelength range from 0.4 to 11.8 mm [
26]. The VIIRS data are from NASA’s ArchiveSet 5200. This article studies the moderate-resolution bands in the visible and NIR region of VIIRS, as shown in
Table 4.
EO-1 was launched in 2000. Hyperion is one of the primary instruments on board and the first space-based high-resolution hyperspectral imaging instrument. It provides 198 calibrated spectral bands, including visible to NIR bands (8–57) and SWIR bands (77–224), covering wavelengths from 430 to 2400 nm. The sampling interval is approximately 10 nm. In this article, Hyperion data were used to provide spectral reflectance profiles [
27].
MODIS images have a large coverage area, high temporal resolution, and the longest operational period among all mentioned satellites, making more angle inputs available for the TOA reflectance model.
Figure 4 shows the band coverage of the four sensors. Differently colored lines represent different sensor bands, with the numbers below indicating their respective band numbers.
Over 7000 Level 1B MODIS images obtained from NASA yield calculated TOA reflectance data and corresponding angle data with which to fit the model parameters. After model establishment, these data were used to validate the accuracy of the long-term TOA reflectance predictions.
To further validate the accuracy and test the computational precision of the TOA reflectance model compared with other satellite sensors, this study used satellite data acquired from EO-1 Hyperion, Landsat 8 OLI, and Suomi NPP VIIRS. OLI’s Level 1T TOA reflectance data and Hyperion’s TOA hyperspectral radiance data are publicly available from the USGS [
28,
29]. The VIIRS data were obtained from NASA’s public database [
30].
This study primarily used Python 3.10 during the data processing and analysis, along with several third-party libraries [
31]. NumPy 1.24 and pandas 2.0 were used for fundamental computations [
32,
33]. Matplotlib 3.7, OpenCV 4.7, and seaborn 0.12 were used for data visualization [
34,
35,
36]. Statsmodels 0.14 was used for modeling [
37]. The acquisition of satellite images was facilitated using the requests 2.31 library. The processing of image data involved the application of pyhdf 0.10, netCDF4 1.6, and GDAL 3.8 libraries [
38,
39]. All of these libraries were adopted for comprehensive data analysis and presentation.
2.3. Data Quality Control
Daily TOA reflectance data were extracted from Terra/Aqua MODIS images between 2003 and 2012 and used to build the TOA reflectance model. A single image may contain cloud cover, reducing the model’s accuracy. Therefore, the data quality was assessed based on the nonuniformity of the target area. Non-clear-sky records were removed by setting a threshold for spatial nonuniformity. Any data with spatial nonuniformity higher than this threshold were considered associated with non-clear-sky weather and were ignored during the model fitting process. The spatial nonuniformity was defined as the coefficient of variation (CV) of all pixels within the target area in the image as follows:
Here,
represents the mean of all pixels within the target area in the satellite image, and
represents the SD. This is a normalized measure of dispersion, also known as the relative standard deviation. An abnormally large CV value indicates cloud cover in the target area, meaning that it is not applicable for radiometric calibration. By calculating the CV for the target area in all images, a histogram of the daily nonuniformity distribution for the Libyan Desert was obtained, as shown in
Figure 5. The CV values roughly followed a normal distribution.
According to the rule, nearly 99.7% of the values drawn from a normal distribution fall within three SDs. Due to the sparsity of the distribution and the presence of outliers, an iterative shrinking method was used to determine the threshold for spatial nonuniformity. This method gradually narrows down the boundary of the CV distribution through iterative steps, approaching the central region of the distribution. In each iteration, and were recalculated, and values outside the interval were removed until all values fell within the interval. In the last iteration, the upper threshold limit was set to . Data points with CV values exceeding this threshold were labeled “non-clear-sky” data and were not used in the modeling process.
The dashed lines in
Figure 5 represent the boundaries of the
intervals for each iteration. For readability, these boundaries are only partially shown. The boundaries for iterations 1–7, which were far from the peak, are not shown, and the boundaries for iterations 9–13, which were densely packed, are not labeled. By the 14th iteration, all values fell within the
interval. At this point, the iteration process was terminated, and the upper limit of 2.4% represented the desired threshold for nonuniformity.
2.4. Sensitivity Analysis of Atmospheric Parameters
During the imaging process of a remote sensor, solar radiation first enters the Earth’s atmosphere from outer space and reaches the Earth’s surface [
40]. It is then reflected from the surface and scattered back to the onboard sensor through the atmosphere. Solar radiation is absorbed and scattered along its transmission path by various atmospheric molecules and aerosols. Due to the influence of atmospheric absorption and scattering, the apparent surface reflectance decreases once the reflected radiation reaches the TOA [
41,
42].
Second Simulation of a Satellite Signal in the Solar Spectrum (6S) is a computer code that accurately simulates the above process [
43,
44]. In this study, the vector version of the 6S model was used to simulate atmospheric radiative transfer in the visible and NIR bands under different atmospheric conditions. By taking geometrical conditions, atmospheric conditions, spectral bands, altitude, aerosol model type, ground reflectance, and other parameters as inputs, this model can produce the apparent reflectance values at the TOA after the traversal of the atmospheric path.
Understanding the influence of the atmosphere is valuable for model development. To quantify the impact of the atmosphere on the radiometric calibration of an onboard sensor, it is necessary to evaluate the sensitivity of the sensor to AOD, water vapor, and ozone. For this purpose, this study considers data from a clear-sky day on 3 July 2019. The atmospheric parameters for that day were extracted from the Aqua MODIS atmospheric product [
24], as shown in
Table 5.
The analysis was performed using a single-variable method. In each test focusing on a specific atmospheric parameter, that parameter was varied while other parameters were held constant. Multiple values were chosen for each parameter based on its distribution. The TOA reflectance was recalculated for each value using the 6S model and compared with the actual TOA reflectance value using
Here,
is the original TOA reflectance calculated by the actual atmospheric parameters (as shown in
Table 5).
represents the reflectance result recalculated after a specific value is taken from the distribution of a single variable among AOD, water vapor, and ozone.
is their relative difference.
After all the calculations were complete, the impacts of different atmospheric parameter variations on radiative transfer were analyzed.
The AOD data over the Libyan Desert were extracted from MODIS MYD08_D3’s Deep_Blue_Aerosol_Optical_Depth_550_Land_Mean dataset. The AOD values were mainly between 0.1 and 0.5, as shown in
Figure 6a.
Nine AOD values ranging from 0.1 to 0.5 at uniform intervals of 0.05 were selected.
Figure 6b shows that in the blue (470 nm) and red (650 nm) bands, the relative difference
varied by 8.0% as the AOD value increased from 0.1 to 0.5. The relative differences for other bands were within 2.3%.
The long-term column water vapor data were extracted from MYD08_D3’s Water_Vapor_Near_Infrared_Clear_Mean dataset. The values over the Libyan Desert ranged mainly between 0.5 and 5 cm, as shown in
Figure 6c.
Ten water vapor values ranging from 0.5 to 5.0 cm were uniformly sampled at intervals of 0.5 cm. As the value increased from 0.5 to 5.0 cm, the relative difference in the 2130 nm SWIR band decreased by 11.6%, while those in the other bands remained below 2.8%, as shown in
Figure 6d.
Considering the electromagnetic absorption of water vapor, it can be inferred that the 2130 nm SWIR band was influenced by atmospheric water vapor absorption. This indicates that the performance of the TOA reflectance model near 2130 nm could also be affected by water vapor [
45].
The total ozone burden data were extracted from the Total_Ozone_Mean dataset in the MYD08_D3 product. The ozone values in the Libyan Desert were primarily distributed between 200 and 350 DU (Dobson unit), as shown in
Figure 6e.
Nine ozone values ranging from 200 to 400 DU were uniformly sampled at intervals of 25 DU.
Figure 6f shows that as the ozone increased from 200 to 400 DU, the relative differences in the green (555 nm) and red (659 nm) bands increased by 5.3% and 4.2%, respectively, indicating the influence of atmospheric ozone absorption [
46]. In other bands, particularly between 865 and 2130 nm, the relative difference remained below 0.01%.
From the single-variable sensitivity analyses above, it is evident that the AOD and ozone had more significant impacts on the visible spectrum, while water vapor significantly influenced the 2130 nm SWIR band. This suggests that atmospheric parameters could affect a model’s accuracy, particularly at solid absorption wavelengths. Therefore, it was necessary to perform atmospheric correction and introduce atmospheric parameters into the TOA reflectance model.