Vicarious Radiometric Calibration of the Hyperspectral Imaging Microsatellites SPARK-01 and-02 over Dunhuang , China

Hao Zhang 1 ID , Bing Zhang 1,2, Zhengchao Chen 1,* and Zhihua Huang 1,3 1 Key Laboratory of Digital Earth Science, Institute of Remote Sensing and Digital Earth, Chinese Academy of Sciences, No. 9 Dengzhuang South Road, Beijing 100094, China; zhanghao612@radi.ac.cn (H.Z.); zb@radi.ac.cn (B.Z.); 15831606243@163.com (Z.H.) 2 University of Chinese Academy of Sciences, No. 19(A) Yuquan Road, Shijingshan District, Beijing 100049, China 3 College of Geoscience and Surveying Engineering, China University of Mining & Technology, No. 11 Xueyuan Road, Haidian District, Beijing 100083, China * Correspondence: chenzc@radi.ac.cn; Tel.: +86-10-8217-8775; Fax: +86-10-8217-8177


Introduction
At 3:22 am UTC on 22 December 2016, two wide-swath pushbroom hyperspectral imaging microsatellites, SPARK-01 and -02, which were manufactured by the Shanghai Engineering Center for Microsatellites, were successfully launched at the Jiuquan satellite launch center by the CZ-2D rocket.The spectrometers on the satellites were developed by the Academy of Opto-electronics, Chinese Academy of Sciences, less than one year previously.SPARK-01 and -02 have spectral ranges of 400-1000 nm, a swath of ~100 km, a spatial resolution of 50 m and 2048 pixels along the cross-track direction.The spectrometers use prisms to split the beam into different bands, and thus, the spectral resolution (or full width at half maximum, FWHM) varies from 1 to 10 nm. Figure 1 shows a schematic of the satellite; major satellite characteristics are described in Table 1.

Introduction
At 3:22 am UTC on 22 December 2016, two wide-swath pushbroom hyperspectral imaging microsatellites, SPARK-01 and -02, which were manufactured by the Shanghai Engineering Center for Microsatellites, were successfully launched at the Jiuquan satellite launch center by the CZ-2D rocket.The spectrometers on the satellites were developed by the Academy of Opto-electronics, Chinese Academy of Sciences, less than one year previously.SPARK-01 and -02 have spectral ranges of 400-1000 nm, a swath of ~100 km, a spatial resolution of 50 m and 2048 pixels along the cross-track direction.The spectrometers use prisms to split the beam into different bands, and thus, the spectral resolution (or full width at half maximum, FWHM) varies from 1 to 10 nm. Figure 1 shows a schematic of the satellite; major satellite characteristics are described in Table 1.Table 1.Main characteristics of SPARK satellite and imaging sensor.

Satellite Characteristic Description
Spectral bands 160 bands ranging from 400 to 1000 nm (151 and 153 valid bands for SPARK-01 and -02, respectively) Swath and spatial resolution 100 km with a resolution of 50 m at nadir viewing Signal-to-noise ratio (SNR) ≥100:1 on average with conditions: solar zenith = 45°, ground reflectance = 0.The SPARK satellites are lightweight and inexpensive.They provide the advantages of fine spectral resolution and large swath.These two hyperspectral satellites can be used for applications such as environmental and disaster monitoring, target detection, and precise classification.They provide basic information to support quantitative applications, resource exploration, and business applications [1].However, due to size, weight, and cost limitations, SPARK-01 and -02 do not have Table 1.Main characteristics of SPARK satellite and imaging sensor.

Satellite Characteristic Description
Spectral bands 160 bands ranging from 400 to 1000 nm (151 and 153 valid bands for SPARK-01 and -02, respectively) Swath and spatial resolution 100 km with a resolution of 50 m at nadir viewing Signal-to-noise ratio (SNR) ≥100:1 on average with conditions: solar zenith = 45 The SPARK satellites are lightweight and inexpensive.They provide the advantages of fine spectral resolution and large swath.These two hyperspectral satellites can be used for applications such as environmental and disaster monitoring, target detection, and precise classification.They provide basic information to support quantitative applications, resource exploration, and business applications [1].However, due to size, weight, and cost limitations, SPARK-01 and -02 do not have on-board calibration systems.Also, complete preflight radiometric calibrations were not performed in the laboratory due to the short manufacture time and prioritization of more urgent tasks before the satellite launch.Only the spectral calibration for each detector was conducted by a monochromator in the laboratory.The spectral response curves followed the Gauss function quite well after data processing, and thus, the central wavelengths and the full-width at half-maximum (FWHM) of the SPARK satellites were determined.The averaged central spectral wavelengths of these two satellites are slightly different (Figure 2).Moreover, the spectral smile effect is minor for SPARK-01 but is evident in SPARK-02 (Figure 3).This aspect should be considered in the data processing flow.Therefore, in-orbit vicarious calibration must be used to transform the satellite data into meaningful physical information.Previous studies used reflectance-, irradiance-, and radiance-based techniques [2,3] to successfully calibrate satellites such as the SPOT HRV [4], Landsat TM/ETM [3,5,6], Airborne Visible and Infrared Spectrometer [7], EO-1 Hyperion [8,9], and FY [10,11], MISR [12], Landsat OLI [13], CBERS-4 [14], and many other optical remote sensors [15].The reflectance-and irradiance-based methods have been compared with cross-calibration methods to derive the calibration coefficients for the BJ-1 microsatellite [16].The results showed the irradiance-based method to be superior to the reflectance-based method, especially under low-visibility atmosphere conditions.In reality, vicarious calibration methods have always been used in combination with the pre-launch calibration and on-board calibration to determine calibration accuracy and monitor the sensor's radiometric stability [17][18][19].Apart from the vicarious calibration methods frequently applied to multispectral remote sensing satellites, some novel methods for hyperspectral sensors have also been proposed in recent years, such as the improved irradiance-based method [20] and supervised vicarious calibration method [21,22].A distinguishing characteristic of a hyperspectral sensor is its high spectral resolution, and spectral smile effect and spectral shift may greatly affect the radiometric accuracy near the atmospheric absorption wavelength regions [7, 23,24].Due to the lack of on-board calibrator and pre-launch radiometric calibration, the in-orbit calibration of the SPARK-01 and -02 satellites was achieved via a calibration experiment performed at the dry Dunhuang site in the Gobi Desert in western China from 28 February to 10 March 2017.In-situ measurements, including both ground reflectance and atmospheric parameters, were also acquired during this calibration period.Two vicarious calibration methods (i.e., reflectance-based and irradiance-based) were used independently to predict the top-of-atmosphere (TOA) radiance (LTOA) using MODTRAN ® 5 software.The vicarious method results were then used to obtain the final SPARK-01 and -02 calibration coefficients.

Calibration Site and Measurements
Three simultaneous measurement datasets from the Dunhuang calibration site were required for the SPARK radiometric calibrations: raw data from the SPARK satellites, surface reflectance measurements, and atmospheric measurements.Also, in order to correct for non-uniform phenomenon detection due to differing detector responses, two more observations were performed around the time when the calibration experiment occurred: a 90° yaw observation or slide slither over the bright desert region during daytime and a dark current observation over the open ocean during nighttime.Use of the 90° yaw observation is efficient for correcting the non-uniform radiometric response among different detectors.This technique has been utilized for Hyperion [9], Quickbird [25], RapidEye [26], and Landsat 8 [27].Using this technique, all the pixels along the cross-track direction would observe nearly the same scene.Owing to the wide swath (~100 km) of the SPARK satellites, it is difficult to find a uniform ground site wider than 100 km to permit normalization of the different responses among pixels in the cross-track direction.Thus, 90° yaw observation is necessary to perform the relative radiometric calibration.The surface reflectance measurements were conducted by a spectroradiometer (FieldSpec-4, ASD Inc., Longmont, CO, USA) one hour before and after the SPARK satellite overpass.The atmospheric measurements were acquired by a CE318 sunphotometer, a Microtops II sunphotometer (Solar Light Company, Inc., Glenside, PA, USA), an irradiance sphere combined with an SVC GER1500 spectrograph and radiosonde balloons.The details are illustrated as follows.

Calibration Site
The Dunhuang calibration site (40°5′32.80″N,94°23′35.78″E) is located on the eastern edge of the Kumutage Penniform Desert, which is in the Gobi Desert in northwestern China, about 35 km west of the city of Dunhuang, Gansu Province.The calibration area is approximately 1.2 km above sea level.The entire vicarious calibration target area (30 km × 30 km) is situated on a stabilized alluvial fan (see Figure 4).The area used for the vicarious calibration measurements for the high-and

Calibration Site and Measurements
Three simultaneous measurement datasets from the Dunhuang calibration site were required for the SPARK radiometric calibrations: raw data from the SPARK satellites, surface reflectance measurements, and atmospheric measurements.Also, in order to correct for non-uniform phenomenon detection due to differing detector responses, two more observations were performed around the time when the calibration experiment occurred: a 90 • yaw observation or slide slither over the bright desert region during daytime and a dark current observation over the open ocean during nighttime.Use of the 90 • yaw observation is efficient for correcting the non-uniform radiometric response among different detectors.This technique has been utilized for Hyperion [9], Quickbird [25], RapidEye [26], and Landsat 8 [27].Using this technique, all the pixels along the cross-track direction would observe nearly the same scene.Owing to the wide swath (~100 km) of the SPARK satellites, it is difficult to find a uniform ground site wider than 100 km to permit normalization of the different responses among pixels in the cross-track direction.Thus, 90 • yaw observation is necessary to perform the relative radiometric calibration.The surface reflectance measurements were conducted by a spectroradiometer (FieldSpec-4, ASD Inc., Longmont, CO, USA) one hour before and after the SPARK satellite overpass.The atmospheric measurements were acquired by a CE318 sunphotometer, a Microtops II sunphotometer (Solar Light Company, Inc., Glenside, PA, USA), an irradiance sphere combined with an SVC GER1500 spectrograph and radiosonde balloons.The details are illustrated as follows.

Calibration Site
The Dunhuang calibration site (40 • 5 32.80"N, 94 • 23 35.78"E) is located on the eastern edge of the Kumutage Penniform Desert, which is in the Gobi Desert in northwestern China, about 35 km west of the city of Dunhuang, Gansu Province.The calibration area is approximately 1.2 km above sea level.The entire vicarious calibration target area (30 km × 30 km) is situated on a stabilized alluvial fan (see Figure 4).The area used for the vicarious calibration measurements for the high-and medium-spatial resolution sensors is approximately 400 m × 400 m and is located in the center of the alluvial fan; the surface is covered by cemented gravels.Several years ago, this calibration site was protected by the addition of protective fens along the edges to form a 500 m × 500 m square region.Figure 4 shows a Landsat/OLI image of the Dunhuang calibration site in which the surrounding fens can just be discerned.The local atmosphere is dry with low aerosol loading, which is beneficial for the calibration experiments.The atmospheric aerosol characteristics at the site are typical of a rural continental location, although some larger particles have been observed, possibly originating from sand dunes located to the northwest [10,11,28,29].Figure 4 shows a Landsat/OLI image of the Dunhuang calibration site in which the surrounding fens can just be discerned.The local atmosphere is dry with low aerosol loading, which is beneficial for the calibration experiments.The atmospheric aerosol characteristics at the site are typical of a rural continental location, although some larger particles have been observed, possibly originating from sand dunes located to the northwest [10,11,28,29].

SPARK Satellite Observations
SPARK-01 and -02 data were acquired over the Dunhuang calibration site at 06:48:30 UTC on 7 March 2017 and at 06:52:32 UTC 28 February 2017, respectively.The dark current data and 90 • yaw data for the relative calibration were acquired on 13 March and 11 March 2017, respectively, for the SPARK-01 satellite and on 27 February and 28 February 2017 for the SPARK-2 satellite.Figure 5 shows the SPARK-01 and -02 raw data; a number of vertical strips are evident.Clouds are evident in the SPARK-2 image over the southern and eastern areas of the calibration site.Although these atmospheric conditions are not ideal for SPARK-02 calibration, the observations over the calibration site were not affected by either clouds or shadows (Figure 5b), and, thus, the calibration results are expected to be comparable.Detailed imaging information for the calibration site and the relative radiometric calibration is listed in Tables 2 and 3, respectively.

SPARK Satellite Observations
SPARK-01 and -02 data were acquired over the Dunhuang calibration site at 06:48:30 UTC on 7 March 2017 and at 06:52:32 UTC 28 February 2017, respectively.The dark current data and 90° yaw data for the relative calibration were acquired on 13 March and 11 March 2017, respectively, for the SPARK-01 satellite and on 27 February and 28 February 2017 for the SPARK-2 satellite.Figure 5 shows the SPARK-01 and -02 raw data; a number of vertical strips are evident.Clouds are evident in the SPARK-2 image over the southern and eastern areas of the calibration site.Although these atmospheric conditions are not ideal for SPARK-02 calibration, the observations over the calibration site were not affected by either clouds or shadows (Figure 5b), and, thus, the calibration results are expected to be comparable.Detailed imaging information for the calibration site and the relative radiometric calibration is listed in Tables 2 and 3, respectively.(a)

Ground Reflectance Measurements
In-situ ground surface reflectance was measured over a 400 × 400 m square region one hour before and after the SPARK satellite overpass.The surface consists of cemented gravels of different colors and sizes (from mm to cm), as well as sand just beneath the gravel (Figure 6a).The measurements were taken by an ASD, Inc. spectroradiometer along a fixed route as shown in Figure 6b.Adjacent measurement points were ~40 m apart, and 10 measurements were taken around each measurement point.As a result, a total of nearly 1000 surface measurements were acquired from the Dunhuang calibration site.Thorough site measurements were repeated several times during the experiment in order to verify the stability of the surface reflectance.Measurements taken under clear atmospheric conditions were examined carefully and, after the exclusion of any erroneous measurements, averaged to produce the average reflectance of the calibration site.The ground reflectance measured on different dates during the experimental period is shown in Figure 7.The ground reflectance is relatively stable on different dates, with differences of less than 2%.The desert reflectance during the experimental period was also measured on 7 March 2017 as shown in Figure 8.These data were used to verify the radiometric calibration results.

Ground Reflectance Measurements
In-situ ground surface reflectance was measured over a 400 × 400 m square region one hour before and after the SPARK satellite overpass.The surface consists of cemented gravels of different colors and sizes (from mm to cm), as well as sand just beneath the gravel (Figure 6a).The measurements were taken by an ASD, Inc. spectroradiometer along a fixed route as shown in Figure 6b.Adjacent measurement points were ~40 m apart, and 10 measurements were taken around each measurement point.As a result, a total of nearly 1000 surface measurements were acquired from the Dunhuang calibration site.Thorough site measurements were repeated several times during the experiment in order to verify the stability of the surface reflectance.Measurements taken under clear atmospheric conditions were examined carefully and, after the exclusion of any erroneous measurements, averaged to produce the average reflectance of the calibration site.The ground reflectance measured on different dates during the experimental period is shown in Figure 7.The ground reflectance is relatively stable on different dates, with differences of less than 2%.The desert reflectance during the experimental period was also measured on 7 March 2017 as shown in Figure 8.These data were used to verify the radiometric calibration results.Remote Sens. 2018, 10, 120 8 of 36

Ground Reflectance Measurements
In-situ ground surface reflectance was measured over a 400 × 400 m square region one hour before and after the SPARK satellite overpass.The surface consists of cemented gravels of different colors and sizes (from mm to cm), as well as sand just beneath the gravel (Figure 6a).The measurements were taken by an ASD, Inc. spectroradiometer along a fixed route as shown in Figure 6b.Adjacent measurement points were ~40 m apart, and 10 measurements were taken around each measurement point.As a result, a total of nearly 1000 surface measurements were acquired from the Dunhuang calibration site.Thorough site measurements were repeated several times during the experiment in order to verify the stability of the surface reflectance.Measurements taken under clear atmospheric conditions were examined carefully and, after the exclusion of any erroneous measurements, averaged to produce the average reflectance of the calibration site.The ground reflectance measured on different dates during the experimental period is shown in Figure 7.The ground reflectance is relatively stable on different dates, with differences of less than 2%.The desert reflectance during the experimental period was also measured on 7 March 2017 as shown in Figure 8.These data were used to verify the radiometric calibration results.

Atmospheric Data
Atmospheric measurements acquired during the experimental period include columnar atmospheric parameters (i.e., AOD and total columnar water vapor, or CWV), the diffuse-to-global irradiance ratio, and the radiosonde vertical profile.A CE318 photometer was used to measure the total AOD and CWV.A total of 5 days of valid data were acquired on 25, 26, and 28 February and 4 and 7 March 2017.The Langley calibration method [30] and a modified calibration method were used for non-water and water absorption channels, respectively, to update the calibration coefficients for the sun measurement channels.These calculations were made with the measurements acquired on 7 March 2017, as the atmosphere was stable and aerosol burden was low.Then, the AOD was calculated in each channel using Beer's law and a spectral response function [31].We used the measured pressure from a barometer and columnar ozone and nitrogen dioxide content from the Ozone Monitoring Instrument (OMI).The CWV was retrieved using a 4-parameter method [32].These parameters, which were retrieved within the SPARK overpass period, were averaged over 15 minutes and used as inputs in the calibration process.The 550 nm channel AOD was calculated via logarithmic interpolation of the 440 nm and 675 nm channel AODs.The 550 nm channel AOD and CWV are shown in Figure 9

Atmospheric Data
Atmospheric measurements acquired during the experimental period include columnar atmospheric parameters (i.e., AOD and total columnar water vapor, or CWV), the diffuse-to-global irradiance ratio, and the radiosonde vertical profile.A CE318 photometer was used to measure the total AOD and CWV.A total of 5 days of valid data were acquired on 25, 26, and 28 February and 4 and 7 March 2017.The Langley calibration method [30] and a modified calibration method were used for non-water and water absorption channels, respectively, to update the calibration coefficients for the sun measurement channels.These calculations were made with the measurements acquired on 7 March 2017, as the atmosphere was stable and aerosol burden was low.Then, the AOD was calculated in each channel using Beer's law and a spectral response function [31].We used the measured pressure from a barometer and columnar ozone and nitrogen dioxide content from the Ozone Monitoring Instrument (OMI).The CWV was retrieved using a 4-parameter method [32].These parameters, which were retrieved within the SPARK overpass period, were averaged over 15 minutes and used as inputs in the calibration process.The 550 nm channel AOD was calculated via logarithmic interpolation of the 440 nm and 675 nm channel AODs.The 550 nm channel AOD and CWV are shown in Figure 9 for the SPARK satellite overpass dates; data influenced by clouds were excluded.Stable atmospheric conditions are indicated by the AOD and water vapor content patterns on 7 March 2017.

Atmospheric Data
Atmospheric measurements acquired during the experimental period include columnar atmospheric parameters (i.e., AOD and total columnar water vapor, or CWV), the diffuse-to-global irradiance ratio, and the radiosonde vertical profile.A CE318 photometer was used to measure the total AOD and CWV.A total of 5 days of valid data were acquired on 25, 26, and 28 February and 4 and 7 March 2017.The Langley calibration method [30] and a modified calibration method were used for non-water and water absorption channels, respectively, to update the calibration coefficients for the sun measurement channels.These calculations were made with the measurements acquired on 7 March 2017, as the atmosphere was stable and aerosol burden was low.Then, the AOD was calculated in each channel using Beer's law and a spectral response function [31].We used the measured pressure from a barometer and columnar ozone and nitrogen dioxide content from the Ozone Monitoring Instrument (OMI).The CWV was retrieved using a 4-parameter method [32].These parameters, which were retrieved within the SPARK overpass period, were averaged over 15 minutes and used as inputs in the calibration process.The 550 nm channel AOD was calculated via logarithmic interpolation of the 440 nm and 675 nm channel AODs.The 550 nm channel AOD and CWV are shown in Figure 9   Measurements from a Microtops II sunphotometer were used to verify the accuracy of the CE318 observations, as shown in Figure 9.The AOD measurements are more accurate, with differences of less than 0.02 between the two instruments.However, the water vapor content differs greatly between the two instruments.We speculate that the calibration coefficients for the Microtops II 940 nm channel need to be updated.A lack of such updates would introduce additional error in water vapor retrievals.Nevertheless, the CE318 results are expected to be more reliable, as the CE318 automatic operation mode is used extensively worldwide.The 550 nm AOD and CWV were averaged over 15 min intervals within the satellite overpasses on 7 March and 28 February 2017; the average AOD and CWV values are 0.1928 and 0.3513 g/cm 2 for 7 March and 0.3476 and 0.5379 g/cm 2 for 28 February.A rural aerosol type was chosen for use in MODTRAN due to the barren Gobi Desert surroundings.Also, the angstrom exponent coefficients derived from the 440 nm to 675 nm channel AOD measurements are 0.75 and 0.3519 for the SPARK-01 and -02 overpass times, respectively.The ozone density was 299 and 305 DU on 7 March and 28 February 2017, respectively; these values were derived from NASA OMI data [33].
In addition, radiosonde balloons were released during SPARK satellite overpass periods to measure the vertical profiles of atmospheric pressure, temperature, and humidity; balloons were released at 05:49:55 and 05:05:08 UTC on 28 February and 7 March 2017, respectively.Figure 10 shows the vertical profiles measured on each date.The variations in pressure and temperature with altitude are similar between the two dates.The humidity changes little at altitudes greater than 5000 m.Measurements from a Microtops II sunphotometer were used to verify the accuracy of the CE318 observations, as shown in Figure 9.The AOD measurements are more accurate, with differences of less than 0.02 between the two instruments.However, the water vapor content differs greatly between the two instruments.We speculate that the calibration coefficients for the Microtops II 940 nm channel need to be updated.A lack of such updates would introduce additional error in water vapor retrievals.Nevertheless, the CE318 results are expected to be more reliable, as the CE318 automatic operation mode is used extensively worldwide.The 550 nm AOD and CWV were averaged over 15 min intervals within the satellite overpasses on 7 March and 28 February 2017; the average AOD and CWV values are 0.1928 and 0.3513 g/cm 2 for 7 March and 0.3476 and 0.5379 g/cm 2 for 28 February.A rural aerosol type was chosen for use in MODTRAN due to the barren Gobi Desert surroundings.Also, the angstrom exponent coefficients derived from the 440 nm to 675 nm channel AOD measurements are 0.75 and 0.3519 for the SPARK-01 and -02 overpass times, respectively.The ozone density was 299 and 305 DU on 7 March and 28 February 2017, respectively; these values were derived from NASA OMI data [33].
In addition, radiosonde balloons were released during SPARK satellite overpass periods to measure the vertical profiles of atmospheric pressure, temperature, and humidity; balloons were released at 05:49:55 and 05:05:08 UTC on 28 February and 7 March 2017, respectively.Figure 10 shows the vertical profiles measured on each date.The variations in pressure and temperature with altitude are similar between the two dates.The humidity changes little at altitudes greater than 5000 m.Measurements from a Microtops II sunphotometer were used to verify the accuracy of the CE318 observations, as shown in Figure 9.The AOD measurements are more accurate, with differences of less than 0.02 between the two instruments.However, the water vapor content differs greatly between the two instruments.We speculate that the calibration coefficients for the Microtops II 940 nm channel need to be updated.A lack of such updates would introduce additional error in water vapor retrievals.Nevertheless, the CE318 results are expected to be more reliable, as the CE318 automatic operation mode is used extensively worldwide.The 550 nm AOD and CWV were averaged over 15 min intervals within the satellite overpasses on 7 March and 28 February 2017; the average AOD and CWV values are 0.1928 and 0.3513 g/cm 2 for 7 March and 0.3476 and 0.5379 g/cm 2 for 28 February.A rural aerosol type was chosen for use in MODTRAN due to the barren Gobi Desert surroundings.Also, the angstrom exponent coefficients derived from the 440 nm to 675 nm channel AOD measurements are 0.75 and 0.3519 for the SPARK-01 and -02 overpass times, respectively.The ozone density was 299 and 305 DU on 7 March and 28 February 2017, respectively; these values were derived from NASA OMI data [33].
In addition, radiosonde balloons were released during SPARK satellite overpass periods to measure the vertical profiles of atmospheric pressure, temperature, and humidity; balloons were released at 05:49:55 and 05:05:08 UTC on 28 February and 7 March 2017, respectively.Figure 10 shows the vertical profiles measured on each date.The variations in pressure and temperature with altitude are similar between the two dates.The humidity changes little at altitudes greater than 5000 m.To acquire the diffuse-to-global irradiance ratio data, an irradiance sphere was used with a SVC GER1500 spectrograph at the calibration site to measure the irradiance at ten minutes intervals throughout the day.Each of the measurements outlined below were repeated three consecutive times.The global solar irradiance (L1) was measured fist, followed by the sky diffuse irradiance (L2), which was assessed with a light barrier.Finally, the global solar irradiance (L3) was determined [20].Figure 11a shows the diffuse-to-global irradiance ratios at 550 nm on the date of the SPARK satellite overpass during the Dunhuang experiment.The smooth diffuse-to-global irradiance ratio curve indicates a very stable atmosphere on 7 March 2017.Figure 11b shows diffuse-to-global irradiance ratio for the entire spectrum at the time of the SPARK satellite overpass; the lower aerosol burden on 7 March 2017 caused lower diffuse-to-global irradiance ratios in comparison to those measured on 28 February 2017.Lastly, the diffuse-to-global irradiance ratios were convolved with the spectral response functions of the corresponding SPARK satellites channels.To acquire the diffuse-to-global irradiance ratio data, an irradiance sphere was used with a SVC GER1500 spectrograph at the calibration site to measure the irradiance at ten minutes intervals throughout the day.Each of the measurements outlined below were repeated three consecutive times.The global solar irradiance (L1) was measured fist, followed by the sky diffuse irradiance (L2), which was assessed with a light barrier.Finally, the global solar irradiance (L3) was determined [20].Figure 11a shows the diffuse-to-global irradiance ratios at 550 nm on the date of the SPARK satellite overpass during the Dunhuang experiment.The smooth diffuse-to-global irradiance ratio curve indicates a very stable atmosphere on 7 March 2017.Figure 11b shows diffuse-to-global irradiance ratio for the entire spectrum at the time of the SPARK satellite overpass; the lower aerosol burden on 7 March 2017 caused lower diffuse-to-global irradiance ratios in comparison to those measured on 28 February 2017.Lastly, the diffuse-to-global irradiance ratios were convolved with the spectral response functions of the corresponding SPARK satellites channels.To acquire the diffuse-to-global irradiance ratio data, an irradiance sphere was used with a SVC GER1500 spectrograph at the calibration site to measure the irradiance at ten minutes intervals throughout the day.Each of the measurements outlined below were repeated three consecutive times.The global solar irradiance (L1) was measured fist, followed by the sky diffuse irradiance (L2), which was assessed with a light barrier.Finally, the global solar irradiance (L3) was determined [20].Figure 11a shows the diffuse-to-global irradiance ratios at 550 nm on the date of the SPARK satellite overpass during the Dunhuang experiment.The smooth diffuse-to-global irradiance ratio curve indicates a very stable atmosphere on 7 March 2017.Figure 11b shows diffuse-to-global irradiance ratio for the entire spectrum at the time of the SPARK satellite overpass; the lower aerosol burden on 7 March 2017 caused lower diffuse-to-global irradiance ratios in comparison to those measured on 28 February 2017.Lastly, the diffuse-to-global irradiance ratios were convolved with the spectral response functions of the corresponding SPARK satellites channels.

Methods
Vertical striping effects were evident in the raw SPARK satellites images and were caused by several factors, such as odd-even detector processing, "smile" effects due to optical aberrations and misalignments, and signal output variations caused by electrical design.In addition, dark stripes caused by bad pixels are evident in the raw SPARK images (Figure 5).Therefore, several pre-processing steps should be conducted before the absolute radiometric calibration, including: (1) bad pixel fixing using the cross-track neighboring two pixels; (2) dark current subtraction; (3) de-smiling by cubic spline interpolation according to the pre-launch spectral calibration; and (4) relative calibration to normalize the different responses among cross-track detectors.The relative calibration procedure was applied to SPARK satellite data to correct different detector responses within a band.The dark current image and 90 • yaw bright image were used to calculate the relative calibration coefficients for each detector.Then, the SPARK radiance over the calibration site was propagated to the top of the atmosphere (TOA) by using the measured ground reflectance and atmospheric parameters.The radiometric calibration coefficients were derived by dividing the predicted TOA radiance from the averaged calibration site digital number (DN) curves.Figure 12 shows the DN curves extracted from SPARK-01 and -02 data after spectral smile correction and relative radiometric calibration.These curves were derived from 6 × 6 pixel averaged values.In the next section, the reflectance-and irradiance-based methods are compared with each other.

Methods
Vertical striping effects were evident in the raw SPARK satellites images and were caused by several factors, such as odd-even detector processing, "smile" effects due to optical aberrations and misalignments, and signal output variations caused by electrical design.In addition, dark stripes caused by bad pixels are evident in the raw SPARK images (Figure 5).Therefore, several pre-processing steps should be conducted before the absolute radiometric calibration, including: (1) bad pixel fixing using the cross-track neighboring two pixels; (2) dark current subtraction; (3) de-smiling by cubic spline interpolation according to the pre-launch spectral calibration; and (4) relative calibration to normalize the different responses among cross-track detectors.The relative calibration procedure was applied to SPARK satellite data to correct different detector responses within a band.The dark current image and 90° yaw bright image were used to calculate the relative calibration coefficients for each detector.Then, the SPARK radiance over the calibration site was propagated to the top of the atmosphere (TOA) by using the measured ground reflectance and atmospheric parameters.The radiometric calibration coefficients were derived by dividing the predicted TOA radiance from the averaged calibration site digital number (DN) curves.Figure 12 shows the DN curves extracted from SPARK-01 and -02 data after spectral smile correction and relative radiometric calibration.These curves were derived from 6 × 6 pixel averaged values.In the next section, the reflectance-and irradiance-based methods are compared with each other.

Relative Radiometric Calibration and Spectral De-Smiling Correction
Enabled by flexible satellite controls and the wide swath (100 km), 90° yaw imaging was performed over the bright desert.This unique imaging method involves turning all detectors to observe nearly the same scene along the orbit direction.The number of rows of the SPARK-01 and -02 90° yaw images exceeds 20,000, and this number is sufficient to normalize the different responses among pixels along the cross-track direction.The average column value for each pixel was used to normalize the differing response behaviors in the given pixel.Figure 13 illustrates the normal and 90° yaw imaging methods.

Relative Radiometric Calibration and Spectral De-Smiling Correction
Enabled by flexible satellite controls and the wide swath (100 km), 90 • yaw imaging was performed over the bright desert.This unique imaging method involves turning all detectors to observe nearly the same scene along the orbit direction.The number of rows of the SPARK-01 and -02 90 • yaw images exceeds 20,000, and this number is sufficient to normalize the different responses among pixels along the cross-track direction.The average column value for each pixel was used to normalize the differing response behaviors in the given pixel.Figure 13 illustrates the normal and 90 • yaw imaging methods.

Dark Current Subtraction
The dark current (DC) image was acquired over the open ocean during nighttime.The signal remaining in the image represents the DC; accurate DC values were calculated from the average along the row direction, which can be expressed as: where i is the column index, j is the row index, k is the band index, DC is the dark current acquired over the ocean during nighttime, DN is the raw data from the SPARK image, DN * is the SPARK data after dark current correction, and B is the dark current averaged along the row direction.

De-Smiling Correction
Cubic spline interpolation was used for de-smiling Hyperion data given its advantage of retaining spectral curve features [24,34].It was also adopted to interpolate the image after dark current subtraction from the central wavelengths for each pixel provided by pre-launch spectral calibration into the average wavelength (Figure 2) of all 2048 pixels.Then, the hyperspectral cube data after dark current subtraction DN* in the original spectral wavelength was transformed into the new cube data DN ** in the average wavelength, and is expressed as follows: where < DN *> and <DN **> represent spectral data at the spatial position (i, j) before and after de-smiling correction, respectively; <λ(j)> is the central wavelength values for j pixel from pre-launch spectral calibration; λ is the average wavelength of all 2048 pixels; and cubic_spline represents the cubic spline interpolation method.

Uniform Normalization
After dark current subtraction, the differing response in each detector was corrected through columnar normalization as follows:

Dark Current Subtraction
The dark current (DC) image was acquired over the open ocean during nighttime.The signal remaining in the image represents the DC; accurate DC values were calculated from the average along the row direction, which can be expressed as: where i is the column index, j is the row index, k is the band index, DC is the dark current acquired over the ocean during nighttime, DN is the raw data from the SPARK image, DN* is the SPARK data after dark current correction, and B is the dark current averaged along the row direction.

De-Smiling Correction
Cubic spline interpolation was used for de-smiling Hyperion data given its advantage of retaining spectral curve features [24,34].It was also adopted to interpolate the image after dark current subtraction from the central wavelengths for each pixel provided by pre-launch spectral calibration into the average wavelength (Figure 2) of all 2048 pixels.Then, the hyperspectral cube data after dark current subtraction DN* in the original spectral wavelength was transformed into the new cube data DN** in the average wavelength, and is expressed as follows: where <DN*> and <DN**> represent spectral data at the spatial position (i, j) before and after de-smiling correction, respectively; <λ(j)> is the central wavelength values for j pixel from pre-launch spectral calibration; λ is the average wavelength of all 2048 pixels; and cubic_spline represents the cubic spline interpolation method.

Uniform Normalization
After dark current subtraction, the differing response in each detector was corrected through columnar normalization as follows: where DN** is the SPARK data averaged along the row direction after dark current correction and de-smiling correction, DN * * * is the average of quantity DN** in all the column and A is the relative radiometric correction coefficient used in the uniform normalization.Theoretically, if the satellite flies with a strict yaw angle of 90 • , the use of the same imaging path for each detector would cause a one-pixel delay between two adjacent detectors.However, the delay distance may be less than one pixel due partially to inexact yaw angle control and partially to minor differences between the ground sample distances (GSDs) along the orbit direction and across the orbit direction.This pattern of delays forms an evident line on the 90 • yaw image (Figure 14a,c); the correction methods are listed in Figure 14b,d.The total delay, in pixels, is easily visually estimated from the 90 • yaw image.This approximate estimation is sufficiently accurate for use because the image row number is quite large, which allows us to ignore minor errors in delay estimations.Equation (3a) can be modified to Equation (4), which applies to a similar imaging path along the orbit direction, through the addition of the delay factor as follows: where N and M represent the total column number and row number for SPARK 90 • yaw image, respectively, D denotes the delay lines, and S is the starting line number to perform the average operation.
Remote Sens. 2018, 10, 120 14 of 36 where DN ** is the SPARK data averaged along the row direction after dark current correction and de-smiling correction, * * * is the average of quantity DN ** in all the column and A is the relative radiometric correction coefficient used in the uniform normalization.
Theoretically, if the satellite flies with a strict yaw angle of 90°, the use of the same imaging path for each detector would cause a one-pixel delay between two adjacent detectors.However, the delay distance may be less than one pixel due partially to inexact yaw angle control and partially to minor differences between the ground sample distances (GSDs) along the orbit direction and across the orbit direction.This pattern of delays forms an evident line on the 90° yaw image (Figure 14a,c); the correction methods are listed in Figure 14b,d.The total delay, in pixels, is easily visually estimated from the 90° yaw image.This approximate estimation is sufficiently accurate for use because the image row number is quite large, which allows us to ignore minor errors in delay estimations.Equation (3a) can be modified to Equation (4), which applies to a similar imaging path along the orbit direction, through the addition of the delay factor as follows: where N and M represent the total column number and row number for SPARK 90° yaw image, respectively, D denotes the delay lines, and S is the starting line number to perform the average operation.

Absolute Radiometric Calibrations
Reflectance-based and irradiance-based methods are widely used for absolute vicarious radiometric calibrations in in-situ experiments.Both methods require accurate measurements of spectral reflectance for the ground target, as well as spectral AOD, vertical columnar water content, and other meteorological parameters.For the reflectance-based method, atmospheric scattering and absorption are computed based on these measurements using MODTRAN 5.In principle, the reflectance-based method aerosol model is assumed based on experience, which may introduce much uncertainty as AOD increases.The irradiance-based method uses the measured data in the reflectance-based method with measurements of the diffuse-to-global spectral irradiance ratio at ground level.This additional measurement helps reduce uncertainty in the aerosol model used for the scattering calculations [35].The principles used to calculate the TOA spectral reflectance in the reflectance-and irradiance-based methods are shown by Equations ( 5) and ( 6), respectively [20].Both methods use Equation ( 7) to transform the TOA spectral reflectance into the TOA radiance.
In Equations ( 5)-(7), θ s is the sun zenith angle, θ v is the view zenith angle of the sensor, and φ v−s is the relative azimuth angle between the view azimuth angle and the sun azimuth angle.ρ t is the measured spectral reflectance of the ground target, and ρ a is the reflectance that corresponds to the atmospheric path radiance (or atmospheric intrinsic reflectance).S is the atmospheric hemisphere reflectance.T(θ s ) and T(θ v ) are the total transmittance of the solar path and the view path, respectively, while ρ * and L are the TOA spectral reflectance and the TOA radiance of the ground target, respectively.µ s and µ v are the values of cos θ s and cos θ v , respectively, and a s and a v are the diffuse-to-global ratios of the sun direction and the view direction, respectively.E 0 is the TOA solar irradiance, and d is the Sun-Earth distance in astronomical units (AU).
If the atmospheric conditions are stable, a linear relationship exists between the relative optical air mass (m) (i.e., inverse of the cosine of the solar zenith (1/µ s )) and the natural logarithm of 1 minus the diffuse-to-global irradiance ratio (ln(1 − a s )) [2].
On stable days, the fitted slope value (1 − b)τ can be used to compute the diffuse-to-global irradiance ratio for both the solar direction and the viewing direction.During the experiment, diffuse-to-global measurements were taken every 10 min throughout the day.Therefore, the α s can be interpolated with sufficient accuracy via the use of an adjacent measurement in Equation (8).However, the α v must be extrapolated to a zenith angle approximating zero from measurements at observation angles quite different from zero.Thus, if the atmosphere was not very stable, as on 28 February 2017, only α s was used to replace the scattering effect in the reflectance-based method; the upward transmittance was also calculated from MODTRAN 5. Similar modifications have been used previously for UAV hyperspectral sensor vicarious calibration [20].
The ratio of ln(1 − α s ) to relative optical air mass (m) (which had values of no more than 6) at 549.89 nm was scattered (see Figure 15) for both measurements taken early in the day on 28 February and 7 March 2017.Two clear outliers measured on 28 February were removed due to the influence of clouds.The measurements on 7 March show a nearly linear relationship, which indicates stable atmospheric conditions.However, non-linear behavior is observed on 28 February.Therefore, the irradiance-based method applied to the 7 March measurements took the form of Equation ( 6), while that applied to the 28 February measurements took the form of Equation (9).For comparison, Equation ( 6) was also applied to the 28 February measurements despite the atmospheric instability.Then, the spectral diffuse-to-global irradiance ratio was convolved with the spectral response functions to derive the band-weighted values; the diffuse-to-global irradiance ratio at the viewing direction was calculated according to Equation (8).Then, linear regression was performed for each band.The goodness-of-fit (R 2 ) values are shown in Figure 16.Linear relationships are evident in each band for the 7 March measurements, but rather lower linear correlations are noted for the 28 February data.The diffuse-to-global irradiance ratios for both SPARK-01 and SPARK-02 are shown in Figure 17 at their calibration site overpass times.
Remote Sens. 2018, 10, 120 16 of 36 The ratio of ln(1 − ) to relative optical air mass (m) (which had values of no more than 6) at 549.89 nm was scattered (see Figure 15) for both measurements taken early in the day on 28 February and 7 March 2017.Two clear outliers measured on 28 February were removed due to the influence of clouds.The measurements on 7 March show a nearly linear relationship, which indicates stable atmospheric conditions.However, non-linear behavior is observed on 28 February.Therefore, the irradiance-based method applied to the 7 March measurements took the form of Equation ( 6), while that applied to the 28 February measurements took the form of Equation (9).For comparison, Equation ( 6) was also applied to the 28 February measurements despite the atmospheric instability.Then, the spectral diffuse-to-global irradiance ratio was convolved with the spectral response functions to derive the band-weighted values; the diffuse-to-global irradiance ratio at the viewing direction was calculated according to Equation (8).Then, linear regression was performed for each band.The goodness-of-fit (R 2 ) values are shown in Figure 16.Linear relationships are evident in each band for the 7 March measurements, but rather lower linear correlations are noted for the 28 February data.The diffuse-to-global irradiance ratios for both SPARK-01 and SPARK-02 are shown in Figure 17 at their calibration site overpass times.The ratio of ln(1 − ) to relative optical air mass (m) (which had values of no more than 6) at 549.89 nm was scattered (see Figure 15) for both measurements taken early in the day on 28 February and 7 March 2017.Two clear outliers measured on 28 February were removed due to the influence of clouds.The measurements on 7 March show a nearly linear relationship, which indicates stable atmospheric conditions.However, non-linear behavior is observed on 28 February.Therefore, the irradiance-based method applied to the 7 March measurements took the form of Equation ( 6), while that applied to the 28 February measurements took the form of Equation (9).For comparison, Equation ( 6) was also applied to the 28 February measurements despite the atmospheric instability.Then, the spectral diffuse-to-global irradiance ratio was convolved with the spectral response functions to derive the band-weighted values; the diffuse-to-global irradiance ratio at the viewing direction was calculated according to Equation (8).Then, linear regression was performed for each band.The goodness-of-fit (R 2 ) values are shown in Figure 16.Linear relationships are evident in each band for the 7 March measurements, but rather lower linear correlations are noted for the 28 February data.The diffuse-to-global irradiance ratios for both SPARK-01 and SPARK-02 are shown in Figure 17 at their calibration site overpass times.

Relative Radiometric Calibration
The relative radiometric calibration coefficients were derived for the SPARK-01 and -02 satellites.The results in terms of blue, green, and red bands are shown in Figures 18 and 19.In total, 32 sub-regions were evident in the SPARK-01 and -02 dark current curves; this number coincides with the design, which features 32 electrical outputs.These coefficients were applied to SPARK images acquired over the calibration site in on 28 February and 7 March 2017.The non-uniformities and variations were largely eliminated after relative radiometric correction using the row-averaged curves (Figure 20).

Relative Radiometric Calibration
The relative radiometric calibration coefficients were derived for the SPARK-01 and -02 satellites.The results in terms of blue, green, and red bands are shown in Figures 18 and 19.In total, 32 sub-regions were evident in the SPARK-01 and -02 dark current curves; this number coincides with the design, which features 32 electrical outputs.These coefficients were applied to SPARK images acquired over the calibration site in on 28 February and 7 March 2017.The non-uniformities and variations were largely eliminated after relative radiometric correction using the row-averaged curves (Figure 20).

Relative Radiometric Calibration
The relative radiometric calibration coefficients were derived for the SPARK-01 and -02 satellites.The results in terms of blue, green, and red bands are shown in Figures 18 and 19.In total, 32 sub-regions were evident in the SPARK-01 and -02 dark current curves; this number coincides with the design, which features 32 electrical outputs.These coefficients were applied to SPARK images acquired over the calibration site in on 28 February and 7 March 2017.The non-uniformities and variations were largely eliminated after relative radiometric correction using the row-averaged curves (Figure 20).

Absolute Radiometric Calibrations
The MODTRAN-simulated radiance calculated using both the reflectance-and irradiance-based methods is shown in Figures 21 and 22, respectively, for the SPARK-01 and -02 satellites.The absolute radiometric calibration is simple to derive by dividing the radiance from the 6 × 6 averaged DN values.The difference between the results from reflectance-and irradiance-based methods does not exceed 6% for the SPARK-01 satellite and shows evident discrepancies in spectral bands <600 nm.However, the differences between the reflectance-and irradiance-based results are greater than 9% for the SPARK-02 satellite in spectral bands <500 nm.These large discrepancies are caused partially by the relatively large AOT (AOT at 550 nm = 0.35) and partially by the unstable weather conditions on 28 February 2017.In comparison, the improved irradiance-based method, which uses only the downward diffuse-to-global irradiance ratio, derived approximately the same radiance for the SPARK-02 on 28 February 2017 as did the reflectance-based method.

Absolute Radiometric Calibrations
The MODTRAN-simulated radiance calculated using both the reflectance-and irradiance-based methods is shown in Figures 21 and 22, respectively, for the SPARK-01 and -02 satellites.The absolute radiometric calibration is simple to derive by dividing the radiance from the 6 × 6 averaged DN values.The difference between the results from reflectance-and irradiance-based methods does not exceed 6% for the SPARK-01 satellite and shows evident discrepancies in spectral bands <600 nm.However, the differences between the reflectance-and irradiance-based results are greater than 9% for the SPARK-02 satellite in spectral bands <500 nm.These large discrepancies are caused partially by the relatively large AOT (AOT at 550 nm = 0.35) and partially by the unstable weather conditions on 28 February 2017.In comparison, the improved irradiance-based method, which uses only the downward diffuse-to-global irradiance ratio, derived approximately the same radiance for the SPARK-02 on 28 February 2017 as did the reflectance-based method.

Absolute Radiometric Calibrations
The MODTRAN-simulated radiance calculated using both the reflectance-and irradiance-based methods is shown in Figures 21 and 22, respectively, for the SPARK-01 and -02 satellites.The absolute radiometric calibration is simple to derive by dividing the radiance from the 6 × 6 averaged DN values.The difference between the results from reflectance-and irradiance-based methods does not exceed 6% for the SPARK-01 satellite and shows evident discrepancies in spectral bands <600 nm.However, the differences between the reflectance-and irradiance-based results are greater than 9% for the SPARK-02 satellite in spectral bands <500 nm.These large discrepancies are caused partially by the relatively large AOT (AOT at 550 nm = 0.35) and partially by the unstable weather conditions on 28 February 2017.In comparison, the improved irradiance-based method, which uses only the downward diffuse-to-global irradiance ratio, derived approximately the same radiance for the SPARK-02 on 28 February 2017 as did the reflectance-based method.

Discussion
Errors associated with in-situ measurements, data processing, and calibration method selection all contribute to uncertainty in satellite calibration coefficients [20].Given that calibration uncertainty sources, such as ground reflectance measurements, inherent code accuracy, etc. have been thoroughly discussed in the literature, we focused our analysis on uncertainties caused by aerosol type assumptions, AOD measurements, water vapor content measurements, atmospheric profile measurements, and satellite image misregistration.Furthermore, the wavelength shift occurring in hyperspectral data would also impose an additional influence on the radiometric calibration accuracy, especially near the atmospheric absorption wavelengths, due to gases like oxygen, water vapor, carbon dioxide, etc.

Uncertainty Due to Aerosol Type Assumptions
The aerosol type used in the Radiative Transfer Model (RTM) introduces great uncertainty in vicarious calibrations, especially in situations in which the AOD is large.It was impractical to measure the vertical distribution of aerosol characteristics during the calibration campaign.However, the actual aerosol type in Dunhuang, which is in arid northwestern China, is generally close to the RURAL and DESERT types described in MODTRAN.To evaluate uncertainty due to aerosol type, three additional aerosol types similar and dissimilar to local conditions (i.e., urban, desert, and maritime) were chosen to replace the rural aerosol type used in the original calculations.The radiance was computed again using these three aerosol types, and the results were compared to those derived using the rural aerosol type.The angstrom exponent coefficients derived from CE318 measurements were also used as inputs for MODTRAN and were 0.75 and 0.3519 for the SPARK-01 and -02 satellites, respectively.The resulting differences in radiance can be used to evaluate uncertainty due to aerosol type (SPARK-01: Figure 23; SPARK-02: Figure 24).The average relative differences are listed in Table 4.

Discussion
Errors associated with in-situ measurements, data processing, and calibration method selection all contribute to uncertainty in satellite calibration coefficients [20].Given that calibration uncertainty sources, such as ground reflectance measurements, inherent code accuracy, etc. have been thoroughly discussed in the literature, we focused our analysis on uncertainties caused by aerosol type assumptions, AOD measurements, water vapor content measurements, atmospheric profile measurements, and satellite image misregistration.Furthermore, the wavelength shift occurring in hyperspectral data would also impose an additional influence on the radiometric calibration accuracy, especially near the atmospheric absorption wavelengths, due to gases like oxygen, water vapor, carbon dioxide, etc.

Uncertainty Due to Aerosol Type Assumptions
The aerosol type used in the Radiative Transfer Model (RTM) introduces great uncertainty in vicarious calibrations, especially in situations in which the AOD is large.It was impractical to measure the vertical distribution of aerosol characteristics during the calibration campaign.However, the actual aerosol type in Dunhuang, which is in arid northwestern China, is generally close to the RURAL and DESERT types described in MODTRAN.To evaluate uncertainty due to aerosol type, three additional aerosol types similar and dissimilar to local conditions (i.e., urban, desert, and maritime) were chosen to replace the rural aerosol type used in the original calculations.The radiance was computed again using these three aerosol types, and the results were compared to those derived using the rural aerosol type.The angstrom exponent coefficients derived from CE318 measurements were also used as inputs for MODTRAN and were 0.75 and 0.3519 for the SPARK-01 and -02 satellites, respectively.The resulting differences in radiance can be used to evaluate uncertainty due to aerosol type (SPARK-01: Figure 23; SPARK-02: Figure 24).The average relative differences are listed in Table 4. Table 4. Average relative differences in radiance for the SPARK-01 and -02 satellites.

Rural and Maritime
Rural and Desert Reflectance-based 6.58% 0.76% 1.46%Both the irradiance-based and improved irradiance-based methods use the diffuse-to-global ratio to minimize the uncertainty associated with aerosol-type assumptions.The average maximum uncertainty of the irradiance-based method used for SPARK-01 on 7 March 2017 is 2%, or less than half of that calculated for the reflectance-based method.The uncertainty in the reflectance-based method due to the aerosol-type assumption increased to 14% for the SPARK-02 satellite on 28 February 2017 because of the relatively large AOD at 550 nm (0.35).The improved irradiance-based method, despite replacing only the downward transmittance features, considerably decreased uncertainty, i.e., by 9%.In reality, the Dunhuang calibration site is surrounded by the Gobi Desert and, thus, the local aerosol is likely fall into the rural or desert types.However, the uncertainty due to aerosol type was conservatively estimated by using the setting of half the difference between the predicted TOA radiance with the urban aerosol type and that with the rural aerosol type.

Uncertainty Due to AOD Measurements
AOD is retrieved from CE318 measurements with a total uncertainty of ~0.01-0.021,which is spectrally dependent and features higher errors in the UV bands [36]; this uncertainty is validated using CE318 and Microtops II measurements in Section 2.4.Therefore, an uncertainty of ±0.02 was added to the 550 nm AOD used for the SPARK-01 and -02 reflectance-based calibrations.For the reflectance-based method, the uncertainty was estimated by comparing the predicted TOA radiance using different AOD values in MODTRAN 5.For the irradiance-based methods (Equations ( 6) and ( 9)), the uncertainty in the predicted TOA radiance can be attributed to both errors in the directly measured transmittance ∕ and errors in the retrieved CE318 measurements.In reality, the transmittance values consist of the CE318 direct retrievals divided by the calibration coefficient for each channel.Thus, the retrieved transmittance uncertainty is a combination of calibration uncertainty from the CE318 calibration coefficient and uncertainty due to the process of interpolating measured transmittance in a few bands into the SPARK satellite bands.The former (calibration) uncertainty is estimated to be ~0.01-0.02(higher in the UV bands) [36], while the latter is 0.5% of the transmittance [37].It is reasonable to set a relative uncertainty of 0.015 for the measured transmittance in the solar direction because the SPARK satellite spectral range spans from the visible to the near-infrared bands, without UV bands.The transmittance uncertainty in the view Both the irradiance-based and improved irradiance-based methods use the diffuse-to-global ratio to the uncertainty associated with aerosol-type assumptions.The average maximum uncertainty of the irradiance-based method used for SPARK-01 on 7 March 2017 is 2%, or less than half of that calculated for the reflectance-based method.The uncertainty in the reflectance-based method due to the aerosol-type assumption increased to 14% for the SPARK-02 satellite on 28 February 2017 because of the relatively large AOD at 550 nm (0.35).The improved irradiance-based method, despite replacing only the downward transmittance features, considerably decreased uncertainty, i.e., by 9%.In reality, the Dunhuang calibration site is surrounded by the Gobi Desert and, thus, the local aerosol is likely fall into the rural or desert types.However, the uncertainty due to aerosol type was conservatively estimated by using the setting of half the difference between the predicted TOA radiance with the urban aerosol type and that with the rural aerosol type.

Uncertainty Due to AOD Measurements
AOD is retrieved from CE318 measurements with a total uncertainty of ~0.01-0.021,which is spectrally dependent and features higher errors in the UV bands [36]; this uncertainty is validated using CE318 and Microtops II measurements in Section 2.4.Therefore, an uncertainty of ±0.02 was added to the 550 nm AOD used for the SPARK-01 and -02 reflectance-based calibrations.For the reflectance-based method, the uncertainty was estimated by comparing the predicted TOA radiance using different AOD values in MODTRAN 5.For the irradiance-based methods (Equations ( 6) and ( 9)), the uncertainty in the predicted TOA radiance can be attributed to both errors in the directly measured transmittance e −τ/µ s and errors in the retrieved CE318 measurements.In reality, the transmittance values consist of the CE318 direct retrievals divided by the calibration coefficient for each channel.Thus, the retrieved transmittance uncertainty is a combination of calibration uncertainty from the CE318 calibration coefficient and uncertainty due to the process of interpolating measured transmittance in a few bands into the SPARK satellite bands.The former (calibration) uncertainty is estimated to be ~0.01-0.02(higher in the UV bands) [36], while the latter is 0.5% of the transmittance [37].It is reasonable to set a relative uncertainty of 0.015 for the measured transmittance in the solar direction because the SPARK satellite spectral range spans from the visible to the near-infrared bands, without UV bands.The transmittance uncertainty in the view direction can be inferred from that in the solar direction by applying the cosine of the view zenith angle.Also, the uncertainty in the retrieved AOD would cause the path radiance and sphere albedo to change with different signs [2].In order to simplify the calculation, the downward and upward direct transmittances, e −τ/µ s and e −τ/µ v , were replaced with T dir (θ s ) and T dir (θ v ).The transmittance and AOD uncertainties were assumed to be independent; thus, the error propagation equations for the TOA reflectance uncertainties using the irradiance-based (Equation ( 6)) and improved irradiance-based (Equation ( 9)) methods can be written as: The uncertainty estimated for SPARK-01 and -02 using the reflectance-and irradiance-based methods is shown in Figure 25.For the reflectance-based method, an AOD uncertainty of 0.02 contributes little (maximum values of 0.6% and 0.7%, respectively, for SPARK-01 and -02) to the total TOA radiance prediction uncertainty.However, the uncertainties for the irradiance-and improved irradiance-based methods appear higher than that for the reflectance-based method.The average and maximum uncertainties are 2.17% and 2.60%, respectively, for SPARK-01 and 1.20% and 1.45% for SPARK-02.The higher uncertainties for the irradiance-and improved irradiance-based methods may be attributed primarily to the direct transmittance uncertainty, which would be partially decreased by the diffuse transmittance uncertainty calculated by MODTRAN, although with the opposite sign.
angle.Also, the uncertainty in the retrieved AOD would cause the path radiance and sphere albedo to change with different signs [2].In order to simplify the calculation, the downward and upward direct transmittances, ∕ and ∕ , were replaced with Tdir( ) and Tdir( ).The transmittance and AOD uncertainties were assumed to be independent; thus, the error propagation equations for the TOA reflectance uncertainties using the irradiance-based (Equation ( 6)) and improved irradiance-based (Equation ( 9)) methods can be written as: ) ) The uncertainty estimated for SPARK-01 and -02 using the reflectance-and irradiance-based methods is shown in Figure 25.For the reflectance-based method, an AOD uncertainty of 0.02 contributes little (maximum values of 0.6% and 0.7%, respectively, for SPARK-01 and -02) to the total TOA radiance prediction uncertainty.However, the uncertainties for the irradiance-and improved irradiance-based methods appear higher than that for the reflectance-based method.The average and maximum uncertainties are 2.17% and 2.60%, respectively, for SPARK-01 and 1.20% and 1.45% for SPARK-02.The higher uncertainties for the irradiance-and improved irradiance-based methods may be attributed primarily to the direct transmittance uncertainty, which would be partially decreased by the diffuse transmittance uncertainty calculated by MODTRAN, although with the opposite sign.

Uncertainty Due to Water Vapor Measurements
The water vapor content retrieval from the CE318 measurements is expected to have an uncertainty of 10%.Therefore, an uncertainty of ±10% was added to the water vapor content retrievals during the SPARK satellite calibration site overpass.The TOA radiance was computed again with MODTRAN, and the difference represents the calibration uncertainty caused by the water vapor measurement (Figure 26).Large uncertainties are apparent in the water vapor absorption bands near 720, 820, and 940 nm.The uncertainties for the water vapor non-absorption bands are lower than 0.2% and, thus, can be omitted.The reflectance-and irradiance-based methods show similar results (Figure 26).The highest values occur in the 940 nm band, amounting to 4.45% and 4.39% for the reflectance-and irradiance-based (or improved irradiance) methods, respectively, in SPARK-01, and 4.17% and 4.04% in SPARK-02.Due to the low water vapor content in arid areas like Dunhuang, the uncertainty caused by the water vapor measurement is relatively small.

Uncertainty Due to Water Vapor Measurements
The water vapor content retrieval from the CE318 measurements is expected to have an uncertainty of 10%.Therefore, an uncertainty of ±10% was added to the water vapor content retrievals during the SPARK satellite calibration site overpass.The TOA radiance was computed again with MODTRAN, and the difference represents the calibration uncertainty caused by the water vapor measurement (Figure 26).Large uncertainties are apparent in the water vapor absorption bands near 720, 820, and 940 nm.The uncertainties for the water vapor non-absorption bands are lower than 0.2% and, thus, can be omitted.The reflectance-and irradiance-based methods show similar results (Figure 26).The highest values occur in the 940 nm band, amounting to 4.45% and 4.39% for the reflectance-and irradiance-based (or improved irradiance) methods, respectively, in SPARK-01, and 4.17% and 4.04% in SPARK-02.Due to the low water vapor content in arid areas like Dunhuang, the uncertainty caused by the water vapor measurement is relatively small.

Uncertainty Due to Atmospheric Profile Measurements
The vertical distributions of temperature, humidity, pressure, and other atmospheric constituents also influence the TOA radiance prediction.In order to explore uncertainty due to the atmospheric profile, the measured radiosonde data used in MODTRAN 5 were replaced with three atmospheric models (i.e., the Mid-Latitude Summer, MS; Mid-Latitude Winter, MW; and 1976 US Standard Atmosphere, US models).The differences in TOA radiance predicted by the three additional atmospheric models and those by the measured radiosonde data represent the uncertainty due to atmospheric profile measurement, as shown in Figure 27.The irradiance-and improved irradiance-based methods show slightly higher uncertainties due to the atmospheric

Uncertainty Due to Atmospheric Profile Measurements
The vertical distributions of temperature, humidity, pressure, and other atmospheric constituents also influence the TOA radiance prediction.In order to explore uncertainty due to the atmospheric profile, the measured radiosonde data used in MODTRAN 5 were replaced with three atmospheric models (i.e., the Mid-Latitude Summer, MS; Mid-Latitude Winter, MW; and 1976 US Standard Atmosphere, US models).The differences in TOA radiance predicted by the three additional atmospheric models and those by the measured radiosonde data represent the uncertainty due to atmospheric profile measurement, as shown in Figure 27.The irradiance-and improved irradiance-based methods show slightly higher uncertainties due to the atmospheric profile than does the reflectance-based method; however, their uncertainties are less than 1.3% in all bands apart from the water vapor absorption bands near 940 nm and 1135 nm.In addition, the MW model appears to be more similar to the radiosonde measurements, as evidenced by the relatively small difference in the radiances predicted using these two inputs.The MS model is likely to represent the actual conditions, considering the location and season of the calibration experiment.Therefore, the maximum differences, which were derived from replacing the radiosonde measurements with US and MW models, were applied in the calibration uncertainty calculations.
profile than does the reflectance-based method; however, their uncertainties are less than 1.3% in all bands apart from the water vapor absorption bands near 940 nm and 1135 nm.In addition, the MW model appears to be more similar to the radiosonde measurements, as evidenced by the relatively small difference in the radiances predicted using these two inputs.The MS model is likely to represent the actual conditions, considering the location and season of the calibration experiment.Therefore, the maximum differences, which were derived from replacing the radiosonde measurements with US and MW models, were applied in the calibration uncertainty calculations."MS," and "MW" refer to the relative difference in predicted radiance derived from replacing radiosonde measurements with these three atmospheric models.

Uncertainty Due to Image Misregistration Errors
To locate the calibration site, the OLI image acquired on 28 February 2017 was used to geo-rectify the SPARK satellite images.The first-order polynomial method was applied with nearest-neighbor resampling around the calibration site to retain the raw DN acquired by the sensors.A one-pixel misregistration around the calibration site is reasonable between the SPARK and OLI images.In addition, the border of the calibration site can be seen in the OLI image due to its 30 m spatial resolution and high radiometric resolution.Therefore, a misregistration of up to two pixels was assumed in the computation of the average DNs at the calibration site.The average DNs and minimum and maximum average DNs determined by shifting the 6 × 6 pixel area by up to two pixels in all directions are shown in Figure 28.The difference between the averaged DNs and the shifted averaged DNs reflect the uncertainty due to image misregistration errors.The differences caused by misregistration are <1.5% in the SPARK-01 495-955 nm spectral range and the SPARK-02 459-995 nm range.The differences are large at the ends of the spectral range due to high noise; data in these ranges are not generally used.Calibration uncertainties caused by the vertical atmospheric profile measurements for SPARK-01 and -02 using the (a) reflectance-and (b) irradiance-based methods, respectively."US", "MS", and "MW" refer to the relative difference in predicted radiance derived from replacing radiosonde measurements with these three atmospheric models.

Uncertainty Due to Image Misregistration Errors
To locate the calibration site, the OLI image acquired on 28 February 2017 was used to geo-rectify the SPARK satellite images.The first-order polynomial method was applied with nearest-neighbor resampling around the calibration site to retain the raw DN acquired by the sensors.A one-pixel misregistration around the calibration site is reasonable between the SPARK and OLI images.In addition, the border of the calibration site can be seen in the OLI image due to its 30 m spatial resolution and high radiometric resolution.Therefore, a misregistration of up to two pixels was assumed in the computation of the average DNs at the calibration site.The average DNs and minimum and maximum average DNs determined by shifting the 6 × 6 pixel area by up to two pixels in all directions are shown in Figure 28.The difference between the averaged DNs and the shifted averaged DNs reflect the uncertainty due to image misregistration errors.The differences caused by misregistration are <1.5% in the SPARK-01 495-955 nm spectral range and the SPARK-02 459-995 nm range.The differences are large at the ends of the spectral range due to high noise; data in these ranges are not generally used.

Uncertainty Due to Spectral Wavelength Shift
Although the central wavelength values were measured for all the 2048 pixels of SPARK-01 and -02 in the laboratory before launch, the wavelength shift may affect the radiometric calibration result in the band near atmospheric gas absorption wavelength.The spectral shifts of Hyperion were estimated to be 0.38-1.39nm at the 760 nm oxygen band by a spectral fitting algorithm compared with the laboratory spectral calibration [8].The same method was also applied to TG-1 hyperspectral imager and the spectral shifts were 2-3 nm, with an uncertainty of 0.3 nm [38].Thus, the spectral fitting algorithm was also used to estimate the spectral shifts for the SPARK satellite.The measured radiance spectrum over the desert was compared with a MODTARN 5-modeled radiance spectrum using the SPARK spectral calibration parameters to derive the spectral shift value.Figure 29a,c show the comparison near the 760 nm oxygen band in the desert (Figure 5) for SPARK-01 and -02, respectively.The MODTRAN 5 spectrum was normalized to match the SPARK-measured radiance level and the spectral wavelength was shifted in 0.1 nm increments.The optimal shifts were estimated to be −0.1 nm for both SPARK-01 and -02.Such a minor spectral shift indicated that SPARK satellites do not undergo an evident spectral shift.Figure 29b,d show the slightly minimized radiance difference after applying a −0.1 nm shift in the SPARK pre-launch laboratory spectral calibration position.

Uncertainty Due to Spectral Wavelength Shift
Although the central wavelength values were measured for all the 2048 pixels of SPARK-01 and -02 in the laboratory before launch, the wavelength shift may affect the radiometric calibration result in the band near atmospheric gas absorption wavelength.The spectral shifts of Hyperion were estimated to be 0.38-1.39nm at the 760 nm oxygen band by a spectral fitting algorithm compared with the laboratory spectral calibration [8].The same method was also applied to TG-1 hyperspectral imager and the spectral shifts were 2-3 nm, with an uncertainty of 0.3 nm [38].Thus, the spectral fitting algorithm was also used to estimate the spectral shifts for the SPARK satellite.The measured radiance spectrum over the desert was compared with a MODTARN 5-modeled radiance spectrum using the SPARK spectral calibration parameters to derive the spectral shift value.Figure 29a,c show the comparison near the 760 nm oxygen band in the desert (Figure 5) for SPARK-01 and -02, respectively.The MODTRAN 5 spectrum was normalized to match the SPARK-measured radiance level and the spectral wavelength was shifted in 0.1 nm increments.The optimal shifts were estimated to be −0.1 nm for both SPARK-01 and -02.Such a minor spectral shift indicated that SPARK satellites do not undergo an evident spectral shift.Figure 29b,d show the slightly minimized radiance difference after applying a −0.1 nm shift in the SPARK pre-launch laboratory spectral calibration position.

Uncertainty Due to Spectral Wavelength Shift
Although the central wavelength values were measured for all the 2048 pixels of SPARK-01 and -02 in the laboratory before launch, the wavelength shift may affect the radiometric calibration result in the band near atmospheric gas absorption wavelength.The spectral shifts of Hyperion were estimated to be 0.38-1.39nm at the 760 nm oxygen band by a spectral fitting algorithm compared with the laboratory spectral calibration [8].The same method was also applied to TG-1 hyperspectral imager and the spectral shifts were 2-3 nm, with an uncertainty of 0.3 nm [38].Thus, the spectral fitting algorithm was also used to estimate the spectral shifts for the SPARK satellite.The measured radiance spectrum over the desert was compared with a MODTARN 5-modeled radiance spectrum using the SPARK spectral calibration parameters to derive the spectral shift value.Figure 29a,c show the comparison near the 760 nm oxygen band in the desert (Figure 5) for SPARK-01 and -02, respectively.The MODTRAN 5 spectrum was normalized to match the SPARK-measured radiance level and the spectral wavelength was shifted in 0.1 nm increments.The optimal shifts were estimated to be −0.1 nm for both SPARK-01 and -02.Such a minor spectral shift indicated that SPARK satellites do not undergo an evident spectral shift.Figure 29b,d Figure 29.Spectral fit result and optimal spectral calibration result with an -0.1 nm shift at the desert in Dunhuang (Figure 5).(a,c) are for SPARK-01 before and after spectral shifting; and (b,d) are for SPARK-02 before and after spectral shifting.
The spectral shift of −0.1 nm was applied to SPARK-01 and -02 to calculate its contribution to radiometric calibration accuracy (Figure 30).As the SPARK satellite uses a prism, the absolute spectral shift for each band will be linear to its FMHW, expressed as: Considering the additional errors caused by the spectral fitting algorithm itself and laboratory calibration, a ±1 nm spectral shift at the 760 nm band was assumed to further estimate the influence on the radiometric calibration of SPARK satellites.The wavelength position at the 760 nm band was shifted by ±1 nm and the radiance difference was calculated for SPARK-01 and -02.As expected, the uncertainty is evident near the atmospheric gas absorption wavelengths, e.g., Fraunhofer 430 nm The spectral shift of −0.1 nm was applied to SPARK-01 and -02 to calculate its contribution to radiometric calibration accuracy (Figure 30).As the SPARK satellite uses a prism, the absolute spectral shift for each band will be linear to its FMHW, expressed as: Considering the additional errors caused by the spectral fitting algorithm itself and laboratory calibration, a ±1 nm spectral shift at the 760 nm band was assumed to further estimate the influence on the radiometric calibration of SPARK satellites.The wavelength position at the 760 nm band was shifted by ±1 nm and the radiance difference was calculated for SPARK-01 and -02.As expected, the uncertainty is evident near the atmospheric gas absorption wavelengths, e.g., Fraunhofer 430 nm and 685 nm, the 760 nm and 690 nm oxygen bands, and the 720 nm, 820 nm, and 940 nm water vapor bands.The maximum value occurred for the 760 nm oxygen band.The uncertainties are less than 2% in the oxygen bands, and less than 1% in the water vapor bands if a spectral shift of −0.1 nm was ignored during SPARK satellite calibration.However, if the spectral shift was increased to 1 nm, the uncertainties would increase considerably in these atmospheric gas absorption bands (e.g., 8% by the spectral shift of +1 nm for SPARK-01 and >10% by the spectral shift of −1 nm for SPARK-02 (Figure 31).In addition, the uncertainty in the 940 nm band for SPARK-02 is higher than that of SPARK-01 due to the larger water vapor content occurring in the daytime for SPARK-02 radiometric calibration (0.35 g/cm 2 for SPARK-01 versus 0.54 g/cm 2 for SPARK-02).
Remote Sens. 2018, 10, 120 29 of 36 and 685 nm, the 760 nm and 690 nm oxygen bands, and the 720 nm, 820 nm, and 940 nm water vapor bands.The maximum value occurred for the 760 nm oxygen band.The uncertainties are less than 2% in the oxygen bands, and less than 1% in the water vapor bands if a spectral shift of −0.1 nm was ignored during SPARK satellite calibration.However, if the spectral shift was increased to 1 nm, the uncertainties would increase considerably in these atmospheric gas absorption bands (e.g., 8% by the spectral shift of +1 nm for SPARK-01 and >10% by the spectral shift of −1 nm for SPARK-02 (Figure 31).In addition, the uncertainty in the 940 nm band for SPARK-02 is higher than that of SPARK-01 due to the larger water vapor content occurring in the daytime for SPARK-02 radiometric calibration (0.35 g/cm 2 for SPARK-01 versus 0.54 g/cm 2 for SPARK-02).and 685 nm, the 760 nm and 690 nm oxygen bands, and the 720 nm, 820 nm, and 940 nm water vapor bands.The maximum value occurred for the 760 nm oxygen band.The uncertainties are less than 2% in the oxygen bands, and less than 1% in the water vapor bands if a spectral shift of −0.1 nm was ignored during SPARK satellite calibration.However, if the spectral shift was increased to 1 nm, the uncertainties would increase considerably in these atmospheric gas absorption bands (e.g., 8% by the spectral shift of +1 nm for SPARK-01 and >10% by the spectral shift of −1 nm for SPARK-02 (Figure 31).In addition, the uncertainty in the 940 nm band for SPARK-02 is higher than that of SPARK-01 due to the larger water vapor content occurring in the daytime for SPARK-02 radiometric calibration (0.35 g/cm 2 for SPARK-01 versus 0.54 g/cm 2 for SPARK-02).

Total Calibration Uncertainty Estimation
Calibration uncertainties caused by other sources are relatively constant and have been discussed previously [20,39,40].The uncertainty due to ground reflectance measurements has been estimated at 2% in field experiments, and this estimation was validated measurements taken on different days during the calibration experiment.The uncertainty in measurement of diffuse-to-global irradiance ratio contributes 2.0% to the total calibration uncertainty [40].The uncertainty due to ozone measurements with an error of 20% was estimated to be 1.3% [39].Thus, because the ozone acquired from OMI has the uncertainty of 4%, it is reasonable to set this uncertainty to 0.6% [41].Although the accuracy of MODTRAN 5 is much improved and comparable to that of the benchmark Line-by-Line Radiative Transfer Model (LBLRTM) [42], this uncertainty is conservatively estimated to be 1%.The uncertainty due to non-Lambertian ground characteristics was estimated at 2% for the Dunhuang calibration site [16].In total, the overall vicarious calibration uncertainty contains uncertainties and errors caused by atmospheric characterization, surface characterization, radiative transfer calculations, and site-average DN calculations [43].The uncertainties discussed above associated with the reflectance-based method are summarized in Table 5 for the SPARK-01 and -02 satellite calibrations; those associated with the irradiance-based method (used for SPARK-01) and the improved irradiance-based method (used for SPARK-02) are summarized in Table 6.The uncertainties associated with different methods are shown for each spectral band of both satellites in Figure 32.Total uncertainty statistics are listed for different spectral ranges in Table 7.For SPARK-01, uncertainties of 4.71 ± 0.34% and 4.11 ± 0.21% were estimated using the reflectance-and irradiance-based methods, respectively.For SPARK-02, uncertainties of 8.12 ± 0.29% and 5.86 ± 0.29% were estimated at >456 nm using the reflectance-and improved irradiance-based methods, respectively.The uncertainty is greatly increased in other spectral ranges due to high image noise.As expected, the uncertainties in both the irradiance-and improved irradiance-based methods are lower than that in the reflectance-based method, especially when the aerosol optical depth is large (e.g., in the SPARK-02 results).However, the irradiance and improved irradiance-based methods depend greatly on the accuracy of the direct transmittance measurements (in Section 5.2) and the diffuse-to-global irradiance ratios; it is therefore important to improve the accuracy of these measurements.processing would not affect most spectral bands to a great degree, but it may introduce artificial features near the atmospheric gas absorption wavelengths when deriving ground surface reflectance through atmospheric correction.Thus, it is strongly recommended that the spectral polishing technique be applied after atmospheric correction to remove the spectral artificial features near the gas absorption wavelengths.
Remote Sens. 2018, 10, 120 32 of 36 Hyperion data [8].The de-smiling technique is always applied to interpolate the raw data from individual spectral positions into the commonly defined spectral central wavelengths.The de-smiling processing would not affect most spectral bands to a great degree, but it may introduce artificial features near the atmospheric gas absorption wavelengths when deriving ground surface reflectance through atmospheric correction.Thus, it is strongly recommended that the spectral polishing technique be applied after atmospheric correction to remove the spectral artificial features near the gas absorption wavelengths.

Preliminary Validation
The calibration coefficients were derived from the irradiance-based method and applied to the SPARK-01 image.Then, the ground reflectance from the 7 March 2017 SPARK-01 image was calculated using measured atmospheric parameters.The retrieved desert reflectance is compared to the in-situ measured reflectance in Figure 8 (Figure 34).Retrieved values are close to the measured values; the discrepancy is within 8% in 500-1000 nm.The difference is partly attributed to the radiometric calibration and partly attributed to the slight terrain fluctuation, inhomogeneous surface and BRDF effect of the desert.Terra MODIS images acquired on 28 February and 7 March 2017 were also used to verify the vicarious calibration methods.The processes used to predict the TOA radiance were similar to those used in the SPARK calibration.Because MODIS has an on-board calibration system, and thus, its calibration accuracy is expected within 3% [44], the MODIS image radiance was taken as a reference to calculate the relative accuracy of the TOA radiance predicted using the vicarious calibration methods (Table 8).For the MODIS image acquired on 28 February, the improved irradiance-based

Preliminary Validation
The calibration coefficients were derived from the irradiance-based method and applied to the SPARK-01 image.Then, the ground reflectance from the 7 March 2017 SPARK-01 image was calculated using measured atmospheric parameters.The retrieved desert reflectance is compared to the in-situ measured reflectance in Figure 8 (Figure 34).Retrieved values are close to the measured values; the discrepancy is within 8% in 500-1000 nm.The difference is partly attributed to the radiometric calibration and partly attributed to the slight terrain fluctuation, inhomogeneous surface and BRDF effect of the desert.
Remote Sens. 2018, 10, 120 32 of 36 Hyperion data [8].The de-smiling technique is always applied to interpolate the raw data from individual spectral positions into the commonly defined spectral central wavelengths.The de-smiling processing would not affect most spectral bands to a great degree, but it may introduce artificial features near the atmospheric gas absorption wavelengths when deriving ground surface reflectance through atmospheric correction.Thus, it is strongly recommended that the spectral polishing technique be applied after atmospheric correction to remove the spectral artificial features near the gas absorption wavelengths.

Preliminary Validation
The calibration coefficients were derived from the irradiance-based method and applied to the SPARK-01 image.Then, the ground reflectance from the 7 March 2017 SPARK-01 image was calculated using measured atmospheric parameters.The retrieved desert reflectance is compared to the in-situ measured reflectance in Figure 8 (Figure 34).Retrieved values are close to the measured values; the discrepancy is within 8% in 500-1000 nm.The difference is partly attributed to the radiometric calibration and partly attributed to the slight terrain fluctuation, inhomogeneous surface and BRDF effect of the desert.Terra MODIS images acquired on 28 February and 7 March 2017 were also used to verify the vicarious calibration methods.The processes used to predict the TOA radiance were similar to those used in the SPARK calibration.Because MODIS has an on-board calibration system, and thus, its calibration accuracy is expected within 3% [44], the MODIS image radiance was taken as a reference to calculate the relative accuracy of the TOA radiance predicted using the vicarious calibration methods (Table 8).For the MODIS image acquired on 28 February, the improved irradiance-based  Terra MODIS images acquired on 28 February and 7 March 2017 were also used to verify the vicarious calibration methods.The processes used to predict the TOA radiance were similar to those used in the SPARK calibration.Because MODIS has an on-board calibration system, and thus, its calibration accuracy is expected within 3% [44], the MODIS image radiance was taken as a reference to calculate the relative accuracy of the TOA radiance predicted using the vicarious calibration methods (Table 8).For the MODIS image acquired on 28 February, the improved irradiance-based method appears superior to the reflectance-based method in the infrared and shortwave infrared spectral bands, but inferior in the third and fourth bands.However, the weather on 28 February was poor, and thus this comparison shows only that the improved irradiance-based method may be appropriate during non-ideal conditions.On 7 March, the atmosphere was stable and aerosol burden was low.Thus, both the reflectance-and irradiance-based methods predicted values approaching those from the MODIS image, with no more than 4% error in the first four bands.In the 1.2 and 1.6 µm bands, the difference between MODIS and the reflectance-based method is larger than that between MODIS and the irradiance-based method, which is likely due to the aerosol-type assumption.The irradiance-based method shows large difference from MODIS in the 2.1 µm bands, which may be attributed to the low signal-to-noise ratio of the in this band.Comparison of TOA radiances predicted by MODIS and the vicarious calibration methods used for the SPARK satellites show that the SPARK calibration methods achieved the accuracy expected.

Conclusions
This study presents the first in-situ vicarious calibration experiments at the Dunhuang site for the SPARK-01 and -02 satellites.Reflectance-, irradiance-, and improved irradiance-based calibration methods were used on images acquired on 7 March and 28 February 2017 by these two satellites.We proposed a 90 • yaw imaging technique for use in the relative calibration method; such methods are very useful for microsatellites without on-board calibration instruments, and especially for satellites with large swaths.An absolute calibration was performed using MODTRAN 5 data, and the methodological and measurement errors in the calibration results were analyzed in detail.Because the SPARK-01 image was acquired during fair weather (e.g., stable atmosphere and low AOD), the calibration uncertainties of the reflectance-and irradiance-based methods are 4.7% and 4.1%, respectively.However, the SPARK-02 image, which was acquired during poor weather, has an uncertainty of 8.12% using the reflectance-based method from 456 to 1000 nm.Under these conditions, the improved irradiance-based method was superior, producing a lower uncertainty of 5.86%.Thus, the additional diffuse-to-global ratio measurements included in the irradianceand improved irradiance-based methods considerably decreases the calibration uncertainty, likely due to its aerosol property assumptions.The improved irradiance-based method is superior to the reflectance-based method under non-ideal atmospheric conditions as it improves the simulated downward transmittance.Although the irradiance-and improved irradiance-based methods are superior to the reflectance-based method on average, the accuracy of the diffuse-to-global ratio measurements may limit the use of these two methods.Indeed, the instrument used to measure the diffuse-to-global ratio has a lower signal-to-noise ratio in the dark blue bands (i.e., <400 nm) and shortwave infrared bands (i.e., >2.1 µm).Moreover, spectral calibration accuracy is a crucial factor to guarantee accurate radiometric calibration.A 1 nm spectral shift for a hyperspectral sensor with a 10 nm spectral resolution would cause as much as a 10% radiometric calibration error near the gas absorption wavelengths.The precise pre-launch spectral calibration in the laboratory as well as the on-orbit monitoring of spectral wavelength shifting are needed.Also, we strongly suggest combining the calibration results derived by the reflectance-and irradiance-(or improved irradiance-) based methods for optimized results.In the future, irradiance-based methods for hyperspectral satellites should be evaluated in more detail by adding spectrally continuous direct transmittance measurements.This could improve calibration accuracy in the gas absorption bands near 940 nm, 1135 nm, 820 nm, etc.

Figure 1 .
Figure 1.Diagram of the SPARK satellite.

Figure 1 .
Figure 1.Diagram of the SPARK satellite.

Figure 4 .
Figure 4. Dunhuang calibration site for medium-high resolution satellites, as illustrated in a Landsat 8/OLI image acquired on 2 February 2017.(a) Scaled subset image.(b) 5× magnification of a portion of the original image, where the red rectangle is the outline of the calibration site.

Figure 4 .
Figure 4. Dunhuang calibration site for medium-high resolution satellites, as illustrated in a Landsat 8/OLI image acquired on 2 February 2017.(a) Scaled subset image; (b) 5× magnification of a portion of the original image, where the red rectangle is the outline of the calibration site.

Figure 5 .
Figure 5. Subsets of SPARK-01 and -02 data acquired over the Dunhuang calibration site featuring pseudo color composited from the 141 (856.60 nm), 111 (648.60 nm) and 84 (550.30nm) bands.(a) SPARK-01 data acquired on 7 March 2017 at 06:48:30 UTC; (b) SPARK-02 data acquired on 28 February 2017 at 06:52:32 UTC.These images were 180 • rotated from the original raw data to maintain the northern and eastern directions on the top and the right hand, respectively.

Figure 6 .
Figure 6.(a) Photo of the calibration site.(b) Schematic of the surface measurement route.

Figure 7 .
Figure 7. Ground reflectance at the Dunhuang calibration site on different dates.Measurements on 7 March and 28 February 2017 correspond to SPARK-01 and SPARK-02, respectively.

Figure 6 .
Figure 6.(a) Photo of the calibration site; (b) Schematic of the surface measurement route.

Figure 6 .
Figure 6.(a) Photo of the calibration site.(b) Schematic of the surface measurement route.

Figure 7 .
Figure 7. Ground reflectance at the Dunhuang calibration site on different dates.Measurements on 7 March and 28 February 2017 correspond to SPARK-01 and SPARK-02, respectively.

Figure 7 .
Figure 7. Ground reflectance at the Dunhuang calibration site on different dates.Measurements on 7 March and 28 February 2017 correspond to SPARK-01 and SPARK-02, respectively.

Figure 8 .
Figure 8. Ground reflectance measured on 7 March 2017 over the desert area south of the Dunhuang calibration site.

Figure 8 .
Figure 8. Ground reflectance measured on 7 March 2017 over the desert area south of the Dunhuang calibration site.

Remote Sens. 2018, 10 , 120 9 of 36 Figure 8 .
Figure 8. Ground reflectance measured on 7 March 2017 over the desert area south of the Dunhuang calibration site.

Figure 11 .
Figure 11.Diffuse-to-global irradiance ratios measured (a) at 550 nm on the date of the SPARK satellite overpass during the Dunhuang experiment, and (b) for the entire spectrum at the time of the SPARK satellite overpass.

Figure 10 .
Figure 10.Vertical profiles of (a) pressure; (b) temperature; and (c) relative humidity measured using radiosondes released on 7 March and 28 February 2017.

Figure 10 .
Figure 10.Vertical profiles of (a) pressure, (b) temperature, and (c) relative humidity measured using radiosondes released on 7 March and 28 February 2017.

Figure 11 .
Figure 11.Diffuse-to-global irradiance ratios measured (a) at 550 nm on the date of the SPARK satellite overpass during the Dunhuang experiment, and (b) for the entire spectrum at the time of the SPARK satellite overpass.

Figure 11 .
Figure 11.Diffuse-to-global irradiance ratios measured (a) at 550 nm on the date of the SPARK satellite overpass during the Dunhuang experiment, and (b) for the entire spectrum at the time of the SPARK satellite overpass.

Figure 13 .
Figure 13.Schematics of (a) normal and (b) 90° yaw imaging methods in the SPARK satellite.

Figure 13 .
Figure 13.Schematics of (a) normal and (b) 90 • yaw imaging methods in the SPARK satellite.

Figure 14 .
Figure 14.90 • yaw images and delay lines in different columns, including the (a) subset image and (b) delay effect correction scheme for SPARK-01 satellite images and the (c) subset image and (d) delay effect correction scheme for SPARK-02 satellite images.

Figure 17 .
Figure 17.Diffuse-to-global irradiance ratio extrapolated and interpolated to the solar direction ( ) and viewing direction ( ) during the satellite overpasses.

Figure 17 .
Figure 17.Diffuse-to-global irradiance ratio extrapolated and interpolated to the solar direction (α s ) and viewing direction (α v ) during the satellite overpasses.

Figure 20 .
Figure 20.Row-averaged values at 650.4 nm, 551.5 nm, and 461.7 nm from images acquired over the calibration site on 7 March and 28 February 2017 using SPARK-01 (a) before and (b) after relative calibration, and at 638.0 nm, 549.5 nm, and 459.0 nm using SPARK-02 (c) before and (d) after relative calibration.

Figure 23 .Figure 23 .
Figure 23.Relative differences in radiance using different aerosol types for the SPARK-01 calibration for 7 March 2017 using (a) the reflectance-based method and (b) the irradiance-based method.

Figure 23 .Figure 24 .Figure 24 .
Figure 23.Relative differences in radiance using different aerosol types for the SPARK-01 calibration for 7 March 2017 using (a) the reflectance-based method and (b) the irradiance-based method.

Figure 24 .
Figure 24.Relative differences in radiance using different aerosol types for the SPARK-02 calibration for 28 February 2017 using (a) the reflectance-based method and (b) the irradiance-based method.

Figure 25 .
Figure 25.Calibration uncertainties caused by the AOD measurements for SPARK-01 and -02 using the (a) reflectance-and (b) irradiance-based (or improved irradiance-based) methods, respectively.

Figure 26 .
Figure 26.Calibration uncertainties caused by the water vapor content measurements for SPARK-01 and -02 using the (a) reflectance-and (b) irradiance-based (or improved irradiance-based) methods, respectively.

Figure 26 .
Figure 26.Calibration uncertainties caused by the water vapor content measurements for SPARK-01 and -02 using the (a) reflectance-and (b) irradiance-based (or improved irradiance-based) methods, respectively.

Figure 27 .
Figure27.Calibration uncertainties caused by the vertical atmospheric profile measurements for SPARK-01 and -02 using the (a) reflectance-and (b) irradiance-based methods, respectively."US," "MS," and "MW" refer to the relative difference in predicted radiance derived from replacing radiosonde measurements with these three atmospheric models.

Figure 27 .
Figure27.Calibration uncertainties caused by the vertical atmospheric profile measurements for SPARK-01 and -02 using the (a) reflectance-and (b) irradiance-based methods, respectively."US", "MS", and "MW" refer to the relative difference in predicted radiance derived from replacing radiosonde measurements with these three atmospheric models.

Figure 28 .
Figure 28.Averaged DNs from 6 × 6 pixel areas in the calibration site and DNs calculated by shifting the 6 × 6 pixel area by up to two pixels in all directions for SPARK-01 and -02.

Figure 28 .
Figure 28.Averaged DNs from 6 × 6 pixel areas in the calibration site and DNs calculated by shifting the 6 × 6 pixel area by up to two pixels in all directions for SPARK-01 and -02.

Remote Sens. 2018, 10 , 120 27 of 36 Figure 28 .
Figure 28.Averaged DNs from 6 × 6 pixel areas in the calibration site and DNs calculated by shifting the 6 × 6 pixel area by up to two pixels in all directions for SPARK-01 and -02.

Figure 29 .
Figure29.Spectral fit result and optimal spectral calibration result with an -0.1 nm shift at the desert in Dunhuang (Figure5).(a,c) are for SPARK-01 before and after spectral shifting; and (b,d) are for SPARK-02 before and after spectral shifting.

Figure 32 .
Figure 32.Total calibration uncertainty estimated for SPARK-01 and -02 using (a) the reflectance-based method and (b) the irradiance-based method.

Figure 32 .
Figure 32.Total calibration uncertainty estimated for SPARK-01 and -02 using (a) the reflectance-based method and (b) the irradiance-based method.

Figure 34 .
Figure 34.Comparison of in-situ desert reflectance measurements with retrieved reflectance from the SPARK-01 image acquired on 7 March 2017.

Figure 34 .
Figure 34.Comparison of in-situ desert reflectance measurements with retrieved reflectance from the SPARK-01 image acquired on 7 March 2017.

Figure 34 .
Figure 34.Comparison of in-situ desert reflectance measurements with retrieved reflectance from the SPARK-01 image acquired on 7 March 2017.

Table 2 .
SPARK image acquisition information at the Dunhuang calibration site.

Table 3 .
SPARK image acquisition information for the relative calibration.

Table 2 .
SPARK image acquisition information at the Dunhuang calibration site.

Table 3 .
SPARK image acquisition information for the relative calibration.

Table 7 .
Average relative differences for various wavelength ranges.

Table 8 .
Differences between the TOA radiance predicted by the vicarious methods (i.e., the Reflectance-, Irradiance-, and Improved irradiance-based methods) and that from MODIS image radiance.