remote sensing Article
Vicarious Radiometric Calibration of the Hyperspectral Imaging Microsatellites SPARK-01 and -02 over Dunhuang, China Hao Zhang 1 1
2 3
*
ID
, Bing Zhang 1,2 , Zhengchao Chen 1, * and Zhihua Huang 1,3
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;
[email protected] (H.Z.);
[email protected] (B.Z.);
[email protected] (Z.H.) University of Chinese Academy of Sciences, No. 19(A) Yuquan Road, Shijingshan District, Beijing 100049, China College of Geoscience and Surveying Engineering, China University of Mining & Technology, No. 11 Xueyuan Road, Haidian District, Beijing 100083, China Correspondence:
[email protected]; Tel.: +86-10-8217-8775; Fax: +86-10-8217-8177
Received: 8 December 2017; Accepted: 12 January 2018; Published: 17 January 2018
Abstract: Two wide-swath hyperspectral imaging microsatellites, SPARK-01 and -02, were launched on 22 December 2016. Radiometric calibration coefficients were determined for these two satellites via a calibration experiment performed from the end of February to the beginning of March 2017 at the high-altitude, homogenous Dunhuang calibration site in the Gobi Desert in China. In-situ measurements, including ground reflectance, direct transmittance, diffuse-to-global irradiance ratio, and radiosonde vertical profile, were acquired. A unique relative calibration procedure was developed using actual satellite images. This procedure included dark current computation and non-uniform correction processes. The former was computed by averaging multiple lines of long strip imagery acquired over open oceans during nighttime, while the latter was computed using images acquired after the adjustment of the satellite yaw angle to 90◦ . This technique was shown to be suitable for large-swath satellite image relative calibration. After relative calibration, reflectance, irradiance, and improved irradiance-based methods were used to conduct absolute radiometric calibrations in order to predict the top-of-atmosphere (TOA) radiance. The SPARK-01 and -02 satellites passed over the calibration site on 7 March and 28 February 2017, during which time fair and non-ideal weather occurred, respectively. Thus, the SPARK-01 calibration coefficient was derived using reflectance- and irradiance-based methods, while that of SPARK -02 was derived using reflectanceand improved irradiance-based methods. The sources of calibration uncertainty, which include aerosol-type assumptions, transmittance measurements, water vapor content retrieval, spectral wavelength shift and satellite image misregistration, were explored in detail for different calibration methods. Using the reflectance and irradiance-based methods, the total uncertainty for SPARK-01 was estimated to be 4.7% and 4.1%, respectively, in the and represent spectral data at the spatial position (i, j) before and after where and represent spectral at thecentral spatial position (i, j) before de-smiling de-smiling correction, respectively; data is the wavelength valuesand forafter j pixel from correction, respectively; is the central wavelength values for j pixel from pre-launch spectral
pre-launch spectral calibration; λ is the average wavelength of all 2048 pixels; and cubic_spline calibration; λ is the average wavelength of all 2048 pixels; and cubic_spline represents the cubic spline represents themethod. cubic spline interpolation method. interpolation 3.1.3. Uniform Normalization After dark current subtraction, the differing response in each detector was corrected through columnar normalization as follows: follows:
DN
***
11 MM DN ∗∗ (i,* j,* k) ( iDN , ∗∗∗ k )(i, k=) = M ∑ D N (i , j , k ) M j =1
(3a) (3a)
j =1
DN
* * *DN
(k
1 NN = 1 ∑ DN ∗∗∗ (*i,*k*) ) = N i =1 D N (i , k )
∗∗∗ ( k )
N
i
=1
(3b) (3b)
Remote Sens. 2018, 10, 120
14 of 36
A(i, k) = DN ∗∗∗ (k)/DN ∗∗∗ (i, k)
(3c)
where DN** is the SPARK data averaged along the row direction after dark current correction and Remote Sens. 2018, 10, 120 14 of 36 de-smiling correction, DN ∗∗∗ is the average of quantity DN** in all the column and A is the relative radiometric correction coefficient used (3c) A(iin ,k the ) = uniform DN ***(k )normalization. / DN ***(i,k ) Theoretically, if the satellite flies with a strict yaw angle of 90◦ , the use of the same imaging where DN ** is the SPARK data averaged along the row direction after dark current correction and path for each detector would cause a one-pixel delay between two adjacent detectors. However, the ∗∗∗ is the average of quantity DN ** in all the column and A is the relative de-smiling correction, delay distance may be less than one due to inexact yaw angle control and partially to radiometric correction coefficientpixel used in thepartially uniform normalization. minor differences between the ground sample distances (GSDs) orbit and across Theoretically, if the satellite flies with a strict yaw angle of 90°,along the usethe of the samedirection imaging path ◦ for each detector would cause a one-pixel delay between two adjacent detectors. However, the delay the orbit direction. This pattern of delays forms an evident line on the 90 yaw image (Figure 14a,c); distance may beare lesslisted than one due14b,d. partially to inexact yaw angle control is and partially to minor the correction methods in pixel Figure The total delay, in pixels, easily visually estimated differences between the ground sample distances (GSDs) along the orbit direction and across ◦ yaw image. This approximate estimation is sufficiently accurate for use becausethe from the 90orbit the image direction. This pattern of delays forms an evident line on the 90° yaw image (Figure 14a,c); the row number is quite large, are which us14b,d. to ignore minor errors in delay estimations. Equation (3a) correction methods listedallows in Figure The total delay, in pixels, is easily visually estimated from theto 90°Equation yaw image. This approximate is sufficiently for usethe because can be modified (4), which appliesestimation to a similar imagingaccurate path along orbitthe direction, image row number is quite large, which allows us to ignore minor errors in delay estimations. through the addition of the delay factor as follows: 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 N − Das follows: 1 DN ∗∗∗ (i, k) = N − ∑ DN ∗∗ (i, j, k) 1 D Nj=− DS(i) * * *** DN (i ,k ) = DN (i ,j ,k )
(4) S(i ) SPARK-01) N −×D i j(=for S(i ) = D/M (4) SS((i )i )==ND−/D/M ( M −SPARK-01) i ) (for SPARK-02) M × i×(for S(i ) = N − D / M × (M − i ) (for SPARK-02) where N and M represent the total column number and row number for SPARK 90◦ yaw where N andD M denotes represent the the total column number forline SPARK 90° yaw image, respectively, delay lines, and Sand is row the number starting number toimage, perform the respectively, D denotes the delay lines, and S is the starting line number to perform the average average operation. operation.
(a)
(b)
(c)
(d)
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.
Remote Sens. 2018, 10, 120
15 of 36
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. ρt × T (θs ) × T (θv ) 1 − ρt × s
(5)
e−τ/µv e−τ/µs × ρ t × (1 − ρ t × s ) × 1 − αs 1 − αv
(6)
ρ∗ (θs , θv , φv−s ) = ρ a (θs , θv , φv−s ) + ρ∗ (θs , θv , φv−s ) = ρ a (θs , θv , φv−s ) +
L = ρ∗ × µs × E0 /(d2 × π )
(7)
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 as and av are the diffuse-to-global ratios of the sun direction and the view direction, respectively. E0 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 − αs ) = ln(1 − ρs) − (1 − b)τm
(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]. ρ∗ (θs , θv , φv−s ) = ρ a (θs , θv , φv−s ) +
ρ × e−τ/µs × T (θv ) 1 − αs
(9)
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
Remote Sens. 2018, 10, 120 Remote Sens. 2018, 10, 120
16 of 36 16 of 36
RemoteThe Sens.ratio 2018, of 10, ln(1 120
− ) to relative optical air mass (m) (which had values of no more than166) of at 36 The ln(1 − (see ) toFigure relative airmeasurements mass (m) (which had values of day no more 6) at 549.89 nmratio was of scattered 15)optical for both taken early in the on 28than February 549.89 nm was scattered (see Figure for bothon measurements takenremoved early in the on 28 February and 7 March 2017. Two clear outliers15) measured 28 February were dueday to the influence of and 7 March 2017. Two clear outliers measured on 28 February were removed due to the influence of influence of clouds. The measurements on 7 March show a nearly linear relationship, which indicates clouds. The measurements on 7 March show a nearly linear relationship, which indicates stable clouds. The measurements on However, 7 March show a nearly linear which indicates stable stable atmospheric conditions. non-linear behavior is relationship, observed 28 February. Therefore, atmospheric conditions. However, non-linear behavior is observed on 28onFebruary. Therefore, the atmospheric conditions. However, non-linear behavior is observed on 28 February. Therefore, the the irradiance-based method applied to the 7 March measurements took the form of Equation (6), irradiance-based method applied to the 7 March measurements took the form of Equation (6), while irradiance-based method applied to the 7 March measurements took the form of Equation (6), while while that applied to the 28 February measurements took the form of Equation (9). For comparison, that applied to the 28 February measurements took the form of Equation (9). For comparison, that applied to the February measurements took the form despite of Equation (9). For comparison, Equation (6) was also applied to the the 28 February February measurements measurements despite the atmospheric atmospheric instability. Equation (6) was also28applied to 28 the instability. Equation (6) was also applied to the 28irradiance February measurements despite the atmospheric Then, diffuse-to-global irradiance ratioratio was convolved with thewith spectral functions Then, the thespectral spectral diffuse-to-global was convolved the response spectralinstability. response Then, the spectral diffuse-to-global irradiance ratio was convolved with the spectral response to derive the band-weighted values; the diffuse-to-global irradiance ratio at the viewing functions to derive the band-weighted values; the diffuse-to-global irradiance ratio at thedirection viewing functions to derive the band-weighted values; the diffuse-to-global irradiance ratio at the viewing was calculated according to Equation (8). Then, linear regression was performed for each direction was calculated according to Equation (8). Then, linear regression was performed forband. each direction calculated according toare Equation (8).Figure Then, linear regression wasevident performed for each The goodness-of-fit (R2 ) values are shown in Figure 16. Linear relationships are in each band. Thewas goodness-of-fit (R2) values shown in 16. Linear relationships are evident inband each band. goodness-of-fit (R2) values are but shown in Figure 16. Linear relationships are each for theThe 7 March but rather lower linear correlations are noted forare the 28evident February data. band for the 7 measurements, March measurements, rather lower linear correlations noted for in the 28 band for the 7 March measurements, but rather lower linear correlations are noted for the 28 The diffuse-to-global irradiance ratios for both SPARK-01 and SPARK-02 are shown in Figure 17 at February data. The diffuse-to-global irradiance ratios for both SPARK-01 and SPARK-02 are shown February data. The diffuse-to-global irradiance ratios for both SPARK-01 and SPARK-02 are shown their calibration site overpass times. in Figure 17 at their calibration site overpass times. in Figure 17 at their calibration site overpass times. 0 0 -0.2 -0.2 -0.4 -0.4 -0.6 -0.6
ln(1-α ln(1-α s) s)
y = -0.2173x + 0.1368 R² = 0.9956 y = -0.2173x + 0.1368 R² = 0.9956
28 Feb. 28 Feb. 7 Mar. 7linear(28 Mar. Feb.) linear(28 Feb.) linear(7 Mar.) linear(7 Mar.)
-0.8 -0.8 -1 -1 0 -1.2 0
1 1
y = -0.199x - 0.0474 y = -0.199x - 0.0474 R² = 0.9865 R² = 0.9865 2 3 2
-1.2
m3 m
4 4
5 5
6 6
7 7
-1.4 -1.4
Figure 15. A scatter plot of ln(1 − ) versus m at 549.89 nm. Figure A scatter scatter plot plot of versus m m at at 549.89 549.89 nm. nm. Figure 15. 15. A of ln(1 ln(1 − − αs)) versus
R2 (%) R2 (%)
100 100 99 99 98 98 97 97 96 96 95 95 94 94 93 93 92 92 91 91 90 90 400 400
SPARK-01 SPARK-01 SPARK-02 SPARK-02
500 500
600 600
700 800 700 800 wavelength (nm) wavelength (nm)
900 900
1000 1000
1100 1100
Figure 16. Goodness-of-fit statistics for diffuse-to-global irradiance ratio measurements according to Figure 16. Goodness-of-fit statistics for diffuse-to-global irradiance ratio measurements according to Equation Figure 16.(8). Goodness-of-fit statistics for diffuse-to-global irradiance ratio measurements according to Equation (8). Equation (8).
17 of 36
Remote Sens. 2018, 10, 120 Remote Sens. 2018, 10, 120
17 of 36 17 of 36
αsαand αvαv s and
Remote Sens. 2018, 10, 120
0.5 0.5 0.45 0.45 0.4 0.4 0.35 0.35 0.3 0.3 0.25 0.25 0.2 0.2 0.15 0.15 0.1 0.1 0.05 0.05 0 0 400 400
for SPARK-01 for SPARK-01 for SPARK-01 for SPARK-01 for SPARK-02 for SPARK-02 for SPARK-02 for SPARK-02
500 500
600 600
700 800 700 800 wavelength (nm) wavelength (nm)
900 900
1000 1000
1100 1100
Figure 17. Diffuse-to-global irradiance ratio extrapolated and interpolated to the solar direction ( ) Figure 17. Diffuse-to-global irradiance irradiance ratio ratio extrapolated extrapolated and andinterpolated interpolatedto tothe thesolar solardirection direction((αs)) Figure 17. Diffuse-to-global and viewing direction ( ) during the satellite overpasses. and viewing direction (α ) during the satellite overpasses. and viewing direction ( v ) during the satellite overpasses.
4. Results 4. Results
coefficient coefficient
4.1. Relative Relative Radiometric Radiometric Calibration Calibration 4.1. 4.1. Relative Radiometric Calibration The relative relative radiometric calibration coefficients were derived for the SPARK-01 and -02 The coefficients werewere derived for thefor SPARK-01 and -02 satellites. The relativeradiometric radiometriccalibration calibration coefficients derived 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, The results termsinofterms blue,ofgreen, and red are shown in Figures 18 18 and total, satellites. Theinresults blue, green, andbands red bands are shown in Figures and19. 19. In In total, 32 sub-regions were evident in the SPARK-01 and -02 dark current curves; this number coincides 32 sub-regions SPARK-01 and -02-02 dark current curves; this number coincides with 32 sub-regions were wereevident evidentininthe the SPARK-01 and dark current curves; this number coincides with the design, which features 32 electrical outputs. These coefficients were applied to SPARK the design, which features 32 electrical outputs. These coefficients were applied to SPARK images 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 acquired over theover calibration site in site on 28 and 7 and March 2017. 2017. The non-uniformities and images acquired the calibration in February on 28 February 7 March The non-uniformities and variations were largely eliminated after relative radiometric correction using the row-averaged variations werewere largely eliminated after relative radiometric correction using the row-averaged curves and variations largely eliminated after relative radiometric correction using the row-averaged curves (Figure 20). (Figure 20). curves (Figure 20). 1.6 1.6 1.4 1.4 1.2 1.2 1 1 0.8 0.8 0.6 0.6 0.4 0.4 0.2 0.2 0 0 0 0
red (#109) red (#109) green (#82) green (#82) blue (#39) blue (#39) 256 256
512 512
768 768
1024
1024 sample sample
(a) (a)
Figure 18. Cont.
1280 1280
1536 1536
1792 1792
2048 2048
Remote Sens. 2018, 10, 120
18 of 36
Remote Sens. 2018, 10, 120
18 of 36
Remote 150 Sens. 2018, 10, 120
DN
DN
145 150 140 145 135 140 130 135 125130 120125 115120 110115 105110 100105 1000
18 of 36
red (#109) green (#82) red (#109) blue (#39) green (#82) blue (#39)
256 0
512
256
512
768 768
1024
sample 1024 sample (b)
1280 1280
1536 1536
1792 1792
2048 2048
(b)
1.6 1.6 1.4 1.4 1.2 1.2 1 1 0.8 0.8 0.60.6
coefficient
coefficient
Figure 18. Relative Figure 18. Relative radiometric radiometric correction correction coefficients coefficients for for SPARK-01 SPARK-01 satellite satellite images images at at 650.4 650.4 nm, nm, Figure 18. Relative radiometric correction coefficients for SPARK-01 satellite images atand 650.4 nm, 551.5 nm, and 461.7 nm showing (a) gain curves (non-uniform correction coefficients) (b) offset 551.5 nm, and 461.7 nm showing (a) gain curves (non-uniform correction coefficients) and (b) offset 551.5 nm, and 461.7 nm showing (a) gain curves (non-uniform correction coefficients) and (b) offset curves curves (dark (dark current). current). curves (dark current).
0.40.4 0.20.2 0 0 0 0
red red(#109) (#109) green green(#84) (#84) blue blue(#40) (#40) 256 256
512 512
768 768
1024 1024
sample sample
1280 1280
1536 1536
1792 1792
2048 2048
1280
1536
1792
2048
(a)
DN DN
150150 145145 140140 135135 130130 125125 120 120 115 115 110 110 105 105 100 100 0 0
red (#109) red (#109) green (#84) green (#84) blue (#40) blue (#40) 256 256
512 512
768 768
1024
1024 sample
sample (b)
(b)
1280
1536
1792
2048
Figure 19. Relative radiometric correction coefficients for SPARK-02 satellite images at 638.0 nm, Figure 19. Relative radiometric radiometric correction coefficients coefficients for for SPARK-02 SPARK-02 satellite satellite images images at at 638.0 638.0 nm, nm, Figure 549.5 19. nm,Relative and 459.0 nm showingcorrection (a) gain curves (non-uniform correction coefficients) and (b) offset 549.5 nm, nm, and and 459.0 459.0 nm nm showing (a) gain curves (non-uniform correction coefficients) and (b) offset 549.5 curves (dark current). showing (a) gain curves (non-uniform correction coefficients) and (b) offset curves (dark current). curves (dark current).
Remote Sens. 2018, 10, 120
19 of 36
Remote Sens. 2018, 10, 120
19 of 36
1400
red (#109) green (#82) blue (#39)
1200
DN
1000 800 600 400 200 0 0
256
512
768
1024
sample
1280
1536
1792
2048
(a) 1400
red (#109) green (#82) blue (#39)
1200
DN
1000 800 600 400 200 0 0
256
512
768
1024
sample
1280
1536
1792
2048
1280
1536
1792
2048
(b) 2400
red (#109) green (#84) blue (#40)
2000
DN
1600 1200 800 400 0 0
256
512
768
1024
sample (c)
Figure 20. Cont.
Remote Sens. 2018, 10, 120
20 of 36
Remote Sens. 2018, 10, 120
20 of 36
2400
Remote Sens. 2018, 10, 120
2000 2400
20 of 36
green (#84) blue (#40)
1600 2000
DN
red (#109) green (#84) blue red(#40) (#109)
1200 1600
DN
800 1200 400800 0400 0 0
256 0
256
512 512
768 768
1024
sample (d)1024
sample
1280 1280
1536 1536
1792 1792
2048 2048
Figure 20. Row-averaged (d)nm, Figure 20. Row-averaged values values at at 650.4 650.4 nm, nm, 551.5 551.5 nm, and and 461.7 461.7 nm nm from from images images acquired acquired over over the the calibration site on 7 March and 28 February 2017 using SPARK-01 (a) before and (b) after relative calibration site on 7 March and 28 February 2017 using SPARK-01 (a) before and (b) after relative Figure 20. Row-averaged values at 650.4 nm, 551.5 nm, and 461.7 nm from images acquired over the calibration, and atat638.0 nm, 549.5 nm,nm, and and 459.0459.0 nm using SPARK-02 (c) before before and (d)and after(d) relative calibration, nm, 549.5 nmSPARK-01 using SPARK-02 after calibrationand site on 638.0 7 March and 28 February 2017 using (a) before(c)and (b) after relative calibration. relative calibration. calibration, and at 638.0 nm, 549.5 nm, and 459.0 nm using SPARK-02 (c) before and (d) after relative calibration.
4.2. Absolute Radiometric Calibrations 4.2. Absolute Radiometric Calibrations 4.2.The Absolute Radiometric Calibrations radiance calculated using both the reflectance- and MODTRAN-simulated The MODTRAN-simulated radiance calculated using both the reflectance- and irradiance-based irradiance-based methods is shown in Figures calculated 21 and 22, respectively, for SPARK-01 and The MODTRAN-simulated radiance using both thethe reflectanceand-02 methods is shown in Figures 21 and 22, respectively, for the SPARK-01 and -02 satellites. The absolute satellites. The absolute radiometric calibration is simple derive by dividing the radiance from the 6 irradiance-based methods is shown in Figures 21 and to 22, for the and -02 radiometric calibration is simple to derive by dividing therespectively, radiance from the SPARK-01 6 × 6 averaged DN × 6satellites. averagedThe DN values.radiometric The difference between the results from irradiance-based absolute calibration is simple to derive byreflectancedividing theand radiance from the 6 values. The difference between the results from reflectance- and irradiance-based methods does × 6 averaged DNexceed values.6% Thefor difference between the results from reflectanceand irradiance-based methods does not the SPARK-01 satellite and shows evident discrepancies in spectral not exceed 6% for the SPARK-01 satellite and shows evident discrepancies in spectral bands
