**3. 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.

**Figure 12.** DN curves from SPARK-01 and -02 at the Dunhuang calibration site, averaged over 6 × 6 pixels.

#### *3.1. 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.

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

#### 3.1.1. 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:

$$B(i,k) = \frac{1}{N} \sum\_{j=1}^{N} DC(i,j,k) \tag{1a}$$

$$DN^\*(i,j,k) = DN(i,j,k) - B(i,k) \tag{1b}$$

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.

## 3.1.2. 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:

$$<\langle DN^{\*\*}(i,j,k) \rangle |\_{i,j} = \text{cubic\\_splite}(\langle DN^\*(i,j,k) \rangle |\_{i,j \prime} \langle \lambda(j) \rangle \langle \sqrt{\lambda} \rangle) \tag{2}$$

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; 6*λ*7 is the average wavelength of all 2048 pixels; and *cubic\_spline* represents the cubic spline interpolation method.

## 3.1.3. Uniform Normalization

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

$$DN^{\*\*\*}(i,k) = \frac{1}{M} \sum\_{j=1}^{M} DN^{\*\*}(i,j,k) \tag{3a}$$

$$\overline{DN^{\ast\ast\ast}}(k) = \frac{1}{N} \sum\_{i=1}^{N} DN^{\ast\ast\ast}(i,k) \tag{3b}$$

$$A(i,k) = \overbrace{DN^{\*\*\*}}^{}(k) / DN^{\*\*\*}(i,k) \tag{3c}$$

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:

$$\begin{aligned} DN^{\ast \ast \ast}(i,k) &= \frac{1}{N-D} \sum\_{j=S(i)}^{N-D} DN^{\ast \ast}(i,j,k) \\ S(i) &= D/M \times i \text{ (for SPARTK-01)} \\ S(i) &= N - D/M \times (M-i) \text{ (for SPARTK-02)} \end{aligned} \tag{4}$$

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.

**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.

#### *3.2. 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.

$$
\rho^\*(\theta\_{\rm s}, \theta\_{\rm v}, \phi\_{\rm v-s}) = \rho\_{\rm d}(\theta\_{\rm s}, \theta\_{\rm v}, \phi\_{\rm v-s}) + \frac{\rho\_{\rm t}}{1 - \rho\_{\rm t} \times s} \times T(\theta\_{\rm s}) \times T(\theta\_{\rm v}) \tag{5}
$$

$$
\rho^\*(\theta\_{\mathfrak{s}}, \theta\_{\mathfrak{v}}, \phi\_{\mathfrak{v}-\mathfrak{s}}) = \rho\_{\mathfrak{s}}(\theta\_{\mathfrak{s}}, \theta\_{\mathfrak{v}}, \phi\_{\mathfrak{v}-\mathfrak{s}}) + \frac{e^{-\tau/\mu\_{\mathfrak{s}}}}{1 - a\_{\mathfrak{s}}} \times \rho\_{\mathfrak{t}} \times (1 - \rho\_{\mathfrak{t}} \times \mathfrak{s}) \times \frac{e^{-\tau/\mu\_{\mathfrak{v}}}}{1 - a\_{\mathfrak{v}}} \tag{6}
$$

$$L = \rho^\* \times \mu\_s \times E\_0 / (d^2 \times \pi) \tag{7}$$

In Equations (5)–(7), *θs* is the sun zenith angle, *θv* is the view zenith angle of the sensor, and *φ<sup>v</sup>*−*<sup>s</sup>* 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 *θ<sup>v</sup>*, respectively, and *as* and *av* 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 − *as*)) [2].

$$
\ln(1 - a\_s) = \ln(1 - \rho s) - (1 - b)\tau m \tag{8}
$$

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].

$$
\rho^\*(\theta\_{s\prime}\theta\_{v\prime}\phi\_{v\prime -s}) = \rho\_a(\theta\_{s\prime}\theta\_{v\prime}\phi\_{v\prime -s}) + \frac{\rho \times e^{-\tau/\mu\_s}}{1 - \alpha\_s} \times T(\theta\_v) \tag{9}
$$

The ratio of ln(1 − *<sup>α</sup>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.

**Figure 15.** A scatter plot of ln(1 − *<sup>α</sup>s*) versus *m* at 549.89 nm.

**Figure 16.** Goodness-of-fit statistics for diffuse-to-global irradiance ratio measurements according to Equation (8).

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