http://orcid.org/0000-0003-0417-4190
Figures
Abstract
During the atmospheric correction of remote sensing data in inland waters, the original Second Simulation of the Satellite Signal in the Solar Spectrum-Vector version (6SV) model does not eliminate the specular reflection of downward skylight radiance at the air-water interface. Thus, we propose a modified version of the 6SV model (M6SV) that does remove reflected skylight at the air-water interface. We apply the new model to the atmospheric correction of a Landsat 8 Operational Land Imager (OLI) image over Taihu Lake, China, where the aerosol optical depth is known. In situ reflectance measurements acquired concurrently with the L8/OLI image are used to validate the performance of the new M6SV algorithm. To further analyze the merits and demerits of M6SV, the model is compared with two short-wave infrared (SWIR)-based atmospheric correction models: the Sea-Viewing Wide Field-of-View Sensor Data Analysis System short-wave infrared (SD-SWIR) model and the Vanhellemont & Ruddick short-wave infrared with a per scene fixed aerosol type (VR-SWIR-F) model. Comparisons of results from all three L8/OLI image atmospheric corrections with the in situ remote sensing reflectance data show that M6SV produces reliable atmospheric corrections in the green and red spectral bands and is an effective alternative for Landsat 8 OLI atmospheric correction in inland waters.
Citation: Lu Z, Li J, Shen Q, Zhang B, Zhang H, Zhang F, et al. (2018) Modification of 6SV to remove skylight reflected at the air-water interface: Application to atmospheric correction of Landsat 8 OLI imagery in inland waters. PLoS ONE 13(8): e0202883. https://doi.org/10.1371/journal.pone.0202883
Editor: Rashid Ansumana, Mercy Hospital, SIERRA LEONE
Received: October 18, 2017; Accepted: August 10, 2018; Published: August 24, 2018
Copyright: © 2018 Lu et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Data Availability: All relevant data are within the paper and its Supporting Information file.
Funding: The work was jointly supported by the National Key Research and Development Program of China (2016YFB0501502), the High Resolution Earth Observation Systems of the National Science and Technology Major Projects (41-Y20A31-9003-15/17 to J.L.), the National Natural Science Foundation of China (No. 91638201 to B.Z., 41571361 to Q.S., 41471308 to J.L.).The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Competing interests: The authors have declared that no competing interests exist.
Introduction
Satellite remote sensing is a cost-effective way to monitor and quantify optical, biological, and ecological processes and phenomena in inland waters at large and transboundary scales. However, signals reaching a sensor over water contain both the desired water-leaving surface features and undesired atmospheric effects caused by absorption and scattering. Thus, atmospheric correction, the manipulation that can remove such undesired effects from sensor received signals, is a crucial procedure for inland waters quality monitoring [1].
Ocean color sensors, including SeaWiFS (Sea-viewing Wide Field-of-view Sensor, 1997–2003), MODIS (Moderate-Resolution Imaging Spectroradiometer, 1999–present), MERIS (MEdium Resolution Imaging Spectrometer, 2002–2012), COCTS (Chinese Ocean Color and Temperature Scanner, 2002–present), and VIIRS (Visible Infrared Imaging Radiometer Suite, 2011–present) [2], usually have high revisit frequency at the expense of reduced spatial resolution (250 to 1200 m). However, the components and optical properties of inland waters are more complex than those of oceans because they vary with location and season. Therefore, imagery from sensors with higher spatial resolution may be important for quality monitoring of inland waters by quantitative remote sensing [3].
Although designed for monitoring land objects, the Landsat satellite series, which consists of eight sensors in operation since 1972 with 30 m spatial resolutions, has been used for more detailed observation or long-term application effectively in coastal [4–6] and inland waters [7–12]. An operational land imager (OLI) was recently launched on Landsat 8 and has been particularly useful for studying inland waters. With improvements in data quality and extensions in spectral coverage, L8/OLI has been readily adopted for aquatic science applications [13–15].
Atmospheric corrections, including those for skylight reflection off of the water surface, are necessary for ocean color remote sensing. While atmospherically contaminated signals contain the path radiance and the desired land- or water-leaving radiance over both land and water, signals over water also include specular reflection of downward skylight radiance off of the air-water interface, sun glitter reflection, and whitecaps. The effects of sun glitter and whitecaps are generally ignored [16], but atmospheric corrections over water should, at least, consider the elimination of specular reflection at the air-water interface.
The 6SV (Second Simulation of the Satellite Signal in the Solar Spectrum-Vector version) atmospheric correction method is based on a physical radiative transfer model (RTM), which features a reliable, specific physical meaning and better generalization. The RTM requires the input of meteorological parameters acquired at the time of the satellite overpass. The source codes for 6SV are freely available and have been widely applied in atmospheric corrections over land. A number of researchers have used 6SV for atmospheric corrections over water [17–19]; several researchers have also compared 6SV to other atmospheric corrections methods over water [6,11]. The application of 6SV in atmospheric corrections over water does involve some challenges; for instance, the atmospheric coefficients calculated by 6SV are intended for surface reflectance (Rs, which is a ratio), not remote sensing reflectance (Rrs, unit: sr-1) [11]. Rrs is the ratio of water-leaving radiance at the air-water interface to downward irradiance and is commonly used in remote sensing over water. The conversion of Rs to Rrs via division by pi is an approximate calculation and does not consider the elimination of specular reflection at the air-water interface.
In this paper, we modify the inland water atmospheric correction algorithm in 6SV to correct for skylight reflected by the water surface. Also, in order to facilitate the retrieval of the nine visible to short-wavelength infrared spectral bands produced by L8/OLI, a new subroutine simulating the L8/OLI measurements is integrated into the algorithm. To examine the performance and applicability of the modified 6SV algorithm, the correction results are compared with those from Sea-Viewing Wide Field-of-View Sensor Data Analysis System short-wave infrared (SD-SWIR) as well as the Vanhellemont & Ruddick short-wave infrared with a per scene fixed aerosol type (VR-SWIR-F). SD-SWIR is the “standard” atmospheric correction algorithm with glint correction [16, 20]. VR-SWIR-F is the modified “standard” atmospheric correction algorithm [21] for extremely turbid waters but without glint correction[3]. For further validation, the corrected results from all three algorithms are compared with in situ reflectance measurements acquired at the time of the satellite overpass.
Materials and methods
2.1 Study area
Taihu Lake is located in eastern China (Fig 1(A)) between 119°53′–120°36′ E and 30°56′–31°33′ N, covers 2338 km2 of Jiangsu Province and Zhejiang Province, and is surrounded by the cities of Wuxi, Huzhou, Yixing, and Suzhou (Fig 1(B)). Taihu Lake has a mean water depth of 1.9 m [12]. The third-largest inland freshwater lake in China, it supplies water to several million residents in nearby cities and plays critical roles in economic development and the regional ecosystem. Taihu Lake drew global attention after a blue-green algae bloom event in 2007, which highlighted the introduction of regional pollution into its waters. According the measured SD (Secchi depth), which measures 0.3 m on average, Taihu Lake can be regarded as an extremely turbid body of water [22].
The geographical locations of (a) Taihu Lake in East China and (b) the in situ measurements samples taken on Oct. 26, 2014 at Taihu Lake. (c) The Rrs spectra acquired from the 13 sites. The figure is for illustrative purposes only.
2.2 In situ Rrs measurements
Because the L8/OLI satellite passed over Taihu Lake at ~2:31 (UTC), the in situ data collection, which covered Gonghu Bay and Meiliang Bay (sites denoted by black dots in Fig 1(B)), was conducted ~2:31 (UTC) over Taihu Lake. The in situ data were collected by research team from the Chinese Academy of Sciences for scientific research only, thus no specific permission was required.
The coordinates of the sampling points were determined using a portable GPS and remote sensing reflectance measurements were taken at each station using a FieldSpec® Pro VNIR, ASD spectroradiometer, which covers a spectral range of 350–1050 nm with a spectral resolution of 3 nm. Reflectance was measured using the above-water method [23]. Ten reference plate, water, and skylight spectra were acquired at each station, and the mean spectrum for each type was calculated within ± 5% error to obtain the characteristic remote sensing reflectance (Rrs). The formula for Rrs is as follows:
(1)
where Ed(0+) is the total surface incident irradiance; Lu (λ) denotes the upwelling radiance measured from the water surface; Lsky (λ) is the skylight radiance; Lp (λ) is the measured reference plate radiance; and ρp (λ) is the reference plate reflectance, which was about 30%. rsky is the surface Fresnel contribution, which was interpolated from the lookup tables created by Mobley from measured angles and wind speeds [24–25]. Fig 1(C) shows the Rrs spectra acquired from Taihu Lake on October 26, 2014.
The aerosol optical depth at 550 nm (AOD550) was measured using a hand-held MICROTOPS II Sun photometer linked to a hand-held GPS; this sun photometer measures direct solar radiation in discrete bands selected from five possible channels (440, 675, 870, 936, and 1020 nm) [26]. Ten sets of AOD measurements were acquired in each band in succession at each station. The β and α (described below) of each set were fitted within ± 3% error using Eq (2) [27]. AOD550 values were then calculated, and the mean value was taken.
(2)
where 0τa(λ) is the AOD at wavelength λ, β is the turbidity coefficient, and α denotes the Ångström exponent.
2.3 L8/OLI data and in situ data matching
The L8/OLI is a push-broom scanner with a swath width of 185 km that covers nine spectral bands. It has eight multispectral bands with 30 m spatial resolution and one pan-chromatic band with 15 m spatial resolution. The central wavelengths of the nine bands are 443, 483, 561, 655, 865, 1609, 2201, 591, and 1373 nm. Compared with the L7/ETM+ (enhanced thematic mapper plus) carried on the Landsat 7 mission, L8/OLI has two additional bands (at 443 and 1373 nm) with narrowed original spectral bands. Because of the longer integration times used in the push-broom scanner, L8/OLI offers SNRs (signal-to-noise ratios) approximately three times higher than those produced by the L7/ETM+ [28]. Furthermore, L8/OLI has better quantization, using 12 instead of 8 bits for radiometric digitization [29].
To evaluate the performance of the M6SV model in removing skylight effects from L8/OLI imagery, we used the model for atmospheric corrections over Taihu Lake. The selected L8/OLI image was acquired at 2:31 (UTC) on October 26, 2014 over Taihu Lake. In order to minimize the effects of temporal and spatial mismatches between satellite and in situ data, the time window was narrowed to ~± 1 h of the Landsat 8 overpass. The in situ data collection, which involved 13 sampling stations (denoted by black dots in Fig 1(B)), was conducted from 01:11 to 03:53 (UTC) over Taihu Lake on October 26, 2014. In order to ensure spatial data consistency, model-measured Rrs spectral data were derived by averaging a 1 × 1 pixel area (with ~0.03 km spatial resolution) surrounding the in situ data location.
2.4 Removing skylight reflectance using the modified 6SV model
The 6SV model is an extended version of 6S that takes into consideration light polarization during the signal transfer process [30]. In order to produce atmospheric corrections, 6SV requires inputs such as the meteorological parameters measured at the time of the satellite overpass; the model outputs the atmospheric correction coefficients xa, xb, and xc [31]. Using these three values, the surface reflectance (Rs) can be calculated as follows:
(3)
where λ is the wavelength and Lt is the top-of-the-atmosphere (TOA) radiance received by a satellite sensor in units of W/(m2·sr·μm).
The originally published 6SV version did not consider the elimination of specular reflection of downward skylight radiance at the air-water interface; 6SV retrieves the surface reflectance (Rs) rather than the remote sensing reflectance (Rrs, which has units of sr-1); these reflectance values are expressed in Eq (4) and Eq (1), respectively:
(4)
Thus, we analyzed the successive orders of scattering (SOS) algorithm with the water signal simulation, modified the multiple scattering calculation process in the 6SV source code, and proposed a modified 6SV (M6SV) model to directly generate outputs including the forward scattering radiance (i.e., the downward sky radiance).
If reflections from sun glitter and whitecaps are omitted, the Lt over water received by a satellite sensor is mainly composed of the upward path radiance (Lp (θv–)), the specular reflection downward skylight radiance on the water surface (Lspec), and the desired tLw as follows:
(5)
where t is the diffuse transmittance of the atmosphere; θv−is the zenith angle for the sensor; and “–” indicates the upward path direction. The 6SV model does not directly calculate Lspec, but does give Lp (θv–), from which Lspec can be calculated as follows [32]:
(6)
where Lspec includes two parts: Lspec_1(λ) = Lp(θv+, λ)r(θv)exp(-τ/cosθv) and Lspec_2(λ) = Lp(θv+, λ)r(θs)exp(-τ/cosθs). In these equations, Lspec_1 is specular reflected skylight, which is diffusely transmitted from solar irradiance; Lspec_2 is the path radiance scattered from specular reflected directly transmitted solar irradiance; τ is the total optical depth; r(θv) and r(θs) is the surface Fresnel contribution, θv and θs are the zenith angles for the sensor and the sun, respectively; and Lp (θv+), which can be calculated by modifying the SOS counting process, refers to forward-scattered radiation from the sun direction θs to the observation direction θv.
In 6SV, the path radiance is calculated through the main function-calling subroutines DISCOM, ATMREF, and OSPOL. The optical thickness, atmospheric scattering, and scattering transmittance are computed by DISCOM and ATMREF. OSPOL is the core of the SOS algorithm; its outputs include the normalized radiation field, which includes the upward path radiance Lp (θv–). A detailed flowchart and parameters for these subroutines can be found in related documentation [31].
The goal was to derive a water signal simulation output from the 6SV main function while maintaining the source code structure and atmospheric correction parameters. Thus, the OSPOL code was modified to output the downward radiance Lp (θv+) at ground level and the upward radiance Lp (θv–) at sensor height. Thus, the new path radiance L*p(θv-, λ) = Lp(θv-, λ) + Lspec(λ) recalculated in the ATMREF code replaces the original path radiance Lp (θv–, λ).
Downward radiance from the bottom of the atmosphere was also added to OSPOL and is calculated by changing the upward radiance from the top of the atmosphere to downward radiance from the bottom of the atmosphere. In the primary scattering radiation calculations, downward and upward radiances for any optical thickness τ are computed as follows [31]:
(7)
where I is the first component of the Stokes vector, which describes the radiation intensity; μ is the cosine of the zenith angle; +μ corresponds to upward radiation; and–μ corresponds to downward radiation, where 1 ≥ μ > 0. The parameters ϕ, ω0, F0, and P are, respectively, the cosine of the azimuth angle, the single scattering albedo, the extraterrestrial solar irradiance, and the scattering phase function.
The definitions of the functions for the optical thickness and upward and downward directions are not changed. To introduce incident light from the bottom of the atmosphere, the downward and upward radiances in the primary scattering radiation term above are modified to:
(8)
where τ1 is the total optical thickness. Furthermore, P conforms to the reciprocal principle via:
(9)
The modifications detailed above were added to 6SV as a new subroutine and named OSPOLTOT. We also modified ATMREF to ATMREFTOT and DISCOM to DISCOMTOT to obtain the M6SV model. ATMREFTOT outputs the Lp (θv+, λ) value computed by the subroutine OSPOLTOT, and DISCOMTOT obtains Lp (θv–, λ) and Lp (θv+, λ) by calling subroutine ATMREFTOT. In M6SV, the new outputs xa′, xb′, and xc′ are atmospheric correction parameters that can be used to retrieve the normalized water-leaving reflectance, nρω, directly. Rrs can then be obtained through Eq (10) [33]:
(10)
2.5 Atmospheric corrections using modified 6SV
While the meaning of the model outputs and results are changed by these modifications, the use of the model remains the same. We ran M6SV with synchronous geometrical image and atmospheric parameters to obtain the atmospheric correction parameters xa′, xb′, and xc′. Then, we substituted these correction parameters into the radiometric corrected image, using Eq (11) to obtain the Rrs of each band.
The L8/OLI spectral response function was not included in the original 6SV. For operational efficiency, we added a subroutine containing the spectral response function and modifying the calling command in the main function in M6SV. The spectral response function is resampled at 2.5 nm intervals in the subroutine. The spectral range of each band spans 0.25 to 4 microns; wavelengths falling outside of the effective spectral range are set to zero.
In addition to the modifications described above, the TOA radiance (Lt) and synchronous image geometrical parameters must be prepared before atmospheric corrections can be performed. Fig 2 shows a flow chart of the atmospheric correction calculations in M6SV. The original L8/OLI L1B images are radiometrically calibrated to Lt using the gain and offset parameters extracted from the image metadata file for each band.
L8/OLI L1B, the original L8/OLI L1B images; Lt, the TOA radiance; Date, the image acquisition date; Time, the image acquisition time; θs, the solar zenith angle; φs, the solar azimuth angle; θ0, the sensor zenith angle; φ0, the sensor azimuth angle; L8/OLI SRF, the L8/OLI spectral response function; AOD550, the aerosol optical depth at 550 nm; xa′, xb′, and xc′, atmospheric correction parameters.
Information about θs (the solar zenith angle), φs (the solar azimuth angle), and the image acquisition date and time can be extracted from the image metadata file. Because L8/OLI views nearly vertically [34–35], θ0 (the sensor zenith angle) and φ0 (the sensor azimuth angle) were set to 0° in this study as per the methods of the United States Geological Survey (USGS) [36]. The mean AOD550 values were measured and calculated during in situ data collection. These input data indicate the most likely atmospheric conditions during image acquisition. The configuration parameters are detailed in Table 1. Using these data, M6SV calculates the atmospheric correction coefficients xa′, xb′, and xc′ of each band separately.
We ran VR-SWIR-F to obtain the Rrs of each band, considering the aerosol type fixed over the study area, using 1609 and 2201 nm for aerosol correction and a threshold on the Rayleigh-corrected reflectance in the 1609 nm for cloud and land masking. To obtain the Rrs of each band using the SD-SWIR model, we considered the aerosol type variable over the study area, chose 1609 and 2201 nm for aerosol correction, performed glint correction and cloud masking, and determined aerosol type per pixel.
2.6 Accuracy assessment
To evaluate the precision of the atmospheric correction, we compare the L8/OLI -derived Rrs values from the three different algorithms with those measured in situ. Synchronous image pixels are selected using the sampling site coordinates. In situ data from measurements carried out within ± 1 h of the L8/OLI overpass are chosen; a total of 13 synchropoints (Fig 1(B)) are used in the model comparison analysis.
The precision evaluation indices used for accuracy assessment include the mean ratio (MR), root mean square error (RMSE), and mean relative error (MRE), which are described by the following equations [37–38]:
(12)
where Rcal,i and Rmea,i refer to the Rrs estimated by the model and measured in situ, respectively;
is the average value of the in situ measurements; and N is the number of samples. MR is the mean ratio value between the model-derived and in situ-measured Rrs for each band, where MR values closer to 1 indicate that the model-derived value is closer to the in situ value and is therefore more accurate.
Results and discussion
3.1 Visual inspection of atmospherically-corrected images
The L8/OLI image of the research area was atmospherically corrected using, alternately, theM6SV, VR-SWIR-F and SD-SWIR models. Rrs images derived by the three models at wavelengths of 483, 561, 655, and 865 nm are shown in Fig 3. These Rrs images show overall Rrs spatial distributions for Taihu Lake.
Results for the atmospheric correction of an L8/OLI image (Taihu Lake, China, October 26, 2014) using the M6SV algorithm, for which (a) Rrs(483), (b) Rrs(561), (c) Rrs(655), and (d) Rrs(865) nm, the VR-SWIR-F algorithm, for which (e) Rrs(483), (f) Rrs(561), (g) Rrs(655), and (h) Rrs(865) nm, and the SD-SWIR algorithm, for which (i) Rrs(483), (j) Rrs(561), (k) Rrs(655), and (l) Rrs(865) nm.
Referring to Fig 3, the values of Rrs produced by the three models grow from 483 nm to 561 nm, and then decrease from 561 nm to 865 nm; this results primarily from the spectrally-dependent contributions of inherent optical properties (IOPs), constituent concentrations, and pure water absorption. The M6SV and SD-SWIR models produce larger Rrs values at each wavelength than does the VR-SWIR-F model. Further, Rrs(561) values over Taihu Lake derived by Chen and Zhang [5, 39–40] are typically greater than 0.03 sr-1; the M6SV-derived Rrs(561) and SD-SWIR-derived Rrs(561) values are consistent with this research (Fig 3B and 3J).
Furthermore, the Rrs images derived by VR-SWIR-F and SD-SWIR are missing sections, as indicated by the white patches in Fig 3(E)–3(L); this is a result of the cloud masking operation [3], which defines Rayleigh-corrected reflectance values in Band 6 (centered at 1609 nm) greater than 0.0215 sr-1 as cloud pixels. However, in this case, the exceedingly high reflectance is caused not by clouds, but by an algal bloom event. Moreover, these sections are not missing from the M6SV results (Fig 3(A)–3(D)). Thus, the VR-SWIR-F and SD-SWIR models are invalid during the algal bloom event, while the M6SV model succeeds in obtaining valid atmospheric correction results.
For visible band, the water optical properties in southern Gonghu Bay show typical clear water characteristics. Whereas Taihu Lake is extremely turbid water, then the bottom radiance contributions to the Rrs values should be taken into consideration while the part of the lake is shallow. From Zhushan Bay to the center of the images, the M6SV Rrs values are low in the 483–655 nm bands but high in the 865 nm band. This phenomenon is caused by the algal bloom, for which the featured bands are 483 nm and 865 nm. The 483 nm band captures the strong absorption from algae, as the value of Rrs(483) decreases with increased Chlorophyll-a during an algal bloom. However, low Rrs(483) values can be also found over clear water due to the combined effects of backscattering and absorption. Thus, the high value of Rrs(865) is also used to identify the bloom, as Rrs(865) increases during bloom events in response to strong scattering by phytoplankton particles.
3.2 Methodological comparison
In Fig 4, the M6SV-derived Rrs, VR-SWIR-F-derived Rrs, and SD-SWIR-derived Rrs are plotted against in situ Rrs measurements at 13 synchropoints. Fig 4 also shows a comparison between in situ Rrs measurements and the mean Rrs values at the 13 observation stations (Fig 1) derived using the three atmospheric correction models. The MR and number of synchropoints that fall within ± 15% error are shown in Table 2. Table 2 also reports the corresponding RMSE and MRE values for each model.
L8/OLI Rrs derived using the (a) M6SV, (b) VR-SWIR-F, and (c) SD-SWIR models plotted against in situ measurements from Taihu Lake at 13 synchropoints. The solid line is a 1:1 line; the dashed lines represent ± 15% error around the 1:1 line. The mean ratio (MR) represents the ratio between the model-derived Rrs and in situ-measured Rrs for each band. (d) A comparison of the mean Rrs spectra from the three models and the in situ measurements at the 13 synchropoints.
Fig 4 and Table 2 indicate that the M6SV model offers improved atmospheric correction performance for the L8/OLI image compared to the SD-SWIR model. The M6SV model under-corrects the image, while the VR-SWIR-F model over-corrects the image.
In Fig 4(A)–4(C), the relationship between the data points and the 1:1 line represents the performance of the corresponding model. In order to contextualize the data, ± 15% error lines are given about the 1:1 line; when the data produced by a specific model fall between these dotted lines, the atmospheric correction error of that model is within ± 15%. The distribution of the M6SV data shows that M6SV outperforms (i.e., features lower error than) SD-SWIR at all bands and VR-SWIR-F at 561 nm, while possessing error similar to VR-SWIR-F at 655 and 865 nm.
Fig 4(D) shows that the Rrs spectra retrieved by all three models fall within the same order of magnitude. However, all of the SD-SWIR Rrs values are larger than the M6SV Rrs values. Additionally, when comparing mean values of Rrs produced by M6SV and VR-SWIR-F (Fig 4(D)), it is clear that M6SV performs well at 561 and 655 nm, while VR-SWIR-F performs well at 483 and 865 nm; the evaluation indices in Table 2 support this conclusion. The MR and the MRE values show that the M6SV model performs best at 561 and 655 nm, while the VR-SWIR-F model performs best at 483 and 865 nm. The number of data points within ± 15% error and the RMSE results show that the M6SV model performances best at 561 nm, while the VR-SWIR-F model performs best at 483 nm; the two methods perform comparably at 655 and 865 nm.
For M6SV, the RMSE decreases from 483 to 865 nm, indicating that the precision of the derived result increases. The MRE decreases from 483 to 561 nm, then increases from 561 to 865 nm, and is less than 22.2% at 561 and 655 nm; these values indicate that the M6SV model performed well in the retrieval of Rrs values for Taihu Lake at 561 and 655 nm. At 865 nm, the reflectance signal in the NIR band is exceedingly small, which may explain why the MRE is greater than 50%.
Typically, in turbid water, signals in the blue band (483 nm) are depressed due to the characteristically high absorption by Chlorophyll-a (Ca) and CDOM in this spectral region. For this reason, the green, red, and NIR bands, rather than the blue band, are used for Ca and Total Suspended Sediment (TSS) retrievals [17, 41–43]. Therefore, results from the M6SV model have practical potential for the quantification of quality in extremely turbid waters.
3.3 Sample comparisons
Three example validation comparisons are shown in Fig 5. The three datasets were obtained at different locations and different times. Fig 5(A)–5(C) show Rrs data collected at the locations indicated in Fig 1(B). The collection times of 01:37, 02:34, and 03:53 (UTC) correspond to approximately 1 hour before, the time of, and 1 hour after the L8/OLI overpass. The L8/OLI-measured Rrs spectral data in Fig 5 were derived by averaging a 1 × 1 pixel area surrounding the in situ data site.
(a-c) Examples of validation comparisons between model-derived Rrs and in situ measurements from Taihu Lake matched sampling points. The latitudes and longitudes of the matched sampling points, times of image and in situ spectra acquisition, and models used to derive Rrs are listed. (d) Values of M6SV-derived Rrs at 4 bands changing with the value of AOD550.
The results in Fig 5 further support the conclusions reached above. First, M6SV shows better atmospheric correction performance than SD-SWIR at all bands. The M6SV model also produces better spectral shapes than does the VR-SWIR-F model. It is also apparent that M6SV performs better than VR-SWIR-F at the time of the L8/OLI overpass. M6SV performance appears to be sensitive to timing, which may be caused by a temporal mismatch between satellite and in situ data. Increased in situ data time resolution may improve the M6SV validation results.
Considering relative stability of other M6SV configuration parameters on time scale, the sensitivity of AOD550 to the M6SV are analyzed and discussed here. The M6SV parameters in Table 1 were kept unchanged except for AOD550 values. The AOD550 values are 0.1, 0.2, 0.3, 0.4, 0.5, 0.58, 0.6, 0.7, 0.8, 0.9, 1.2 and 1.5, respectively, for atmospheric correction. For sample (31.37° N, 120.21° E), Rrs values derived by the M6SV at wavelengths of 483, 561, 655, and 865 nm, changing with the AOD550 values, are shown in Fig 5(D). The influence trend in atmospheric correction results of M6SV is basically consistent. That is, the Rrs values decrease linearly with the increase of AOD550 for AOD550 less than 0.6. When AOD550 greater than 0.6, the rate of decrease in Rrs value is increasing with the increase of AOD550. The influence of AOD550 on Rrs decreases with the increase of wavelength. The blue band is the most affected while the near-infrared band is largely unaffected. It may be because attenuation of scattering appears when the wavelength is greater than the aerosol diameter.
3.4 Transect analysis
Rrs spatial distributions derived by M6SV, VR-SWIR-F and SD-SWIR are exhibited in Fig 3. To further demonstrate the differences among the three models, the Rrs values of pixels lying along the transect ‘T’ (Fig 3(E)) are extracted at wavelengths of 483, 561, 655, and 865 nm. These profiles are shown in Fig 6, where the different colors represent the three models used. From left to right, the x-axis pixel values correspond to Gonghu Bay, then Meilang Bay, and finally Zhushan Bay; the selected transect covers turbid waters and the algal bloom area.
Rrs (λ) were derived using the (green) M6SV, (yellow) VR-SWIR-F and (blue) SD-SWIR models at wavelengths of (a) 483, (b) 561, (c) 655, and (d) 865 nm.
In Fig 6, the spectral shapes produced by the three models are generally consistent. The Rrs values derived by SD-SWIR are larger than those derived by M6SV while the Rrs values derived by M6SV are larger than those derived by VR-SWIR-F. With the exception of pixels over the algal bloom, the degree of turbidity along the transect is similar because the Rrs profiles are similar. The 483 and 865 nm bands are most useful for the algal bloom pixels, as explained previously. Thus, for pixels over the algal bloom, the Rrs values derived by M6SV are small at 483 nm and large at 865 nm, while the Rrs values derived by VR-SWIR-F and SD-SWIR are NaNs (Not-a-Number), as these pixels were masked by the algorithms. For pixels over low algal bloom concentrations, the pixels are not masked by the VR-SWIR-F and SD-SWIR model; for such pixels, the Rrs values derived by the two models are small at 483 nm and slightly larger at 865 nm.
3.5 Density-sliced scatterplots
In order to compare the performance of the two models for each pixel over the entire lake area, a L8/OLI image taken over Taihu Lake on October 26, 2014 was processed using both the M6SV and VR-SWIR-F models. Density-sliced scatterplots of the M6SV estimates versus the VR-SWIR-F estimates at wavelengths of 483, 561, 655, and 865 nm are shown in Fig 7. R2 values for the two models range from 0.814 at 483 nm to 0.938 at 655 nm, indicating high correlations. The Rrs values derived by M6SV are larger than those derived by VR-SWIR-F for each band over the entire lake area; this comparison of Rrs values and distribution trends supports the validation and transect results.
Rrs (λ) density-sliced scatterplots at wavelengths of (a) 483, (b) 561, (c) 655, and (d) 865 nm. The plots show the relationship between the M6SV and VR-SWIR-F models for a L8/OLI image taken over Taihu Lake, China, on October 26, 2014.
Conclusions
In this study, a modified 6SV atmospheric correction algorithm (M6SV) is developed to implement a skylight correction. This algorithm is then applied to the atmospheric correction of a L8/OLI image over the extremely turbid inland waters of Taihu Lake in China. To validate the algorithm, M6SV-derived, SD-SWIR-derived, and VR-SWIR-F-derived Rrs values are calculated and compared with the in situ measured reflectance. These comparisons show that the M6SV method improves L8/OLI atmospheric correction performance over the SD-SWIR method at all bands and over the VR-SWIR-F method at 561 and 655 nm. The M6SV MRE values at 561 and 655 nm are less than 22.2%, indicating that the model works well. With the derived Rrs, we can retrieve water quality parameters (such as Ca and TSS) that can be used to monitor the optical and biological properties of inland waters and establish a real-time, wide-ranging, and long-term water quality time series.
Supporting information
S1 Table. Full dataset for the in situ measurements, along with geographical coordinates for the sampling stations.
https://doi.org/10.1371/journal.pone.0202883.s001
(XLSX)
Acknowledgments
We thank USGS for the Landsat 8 OLI imagery and the staff of the Institute of Remote Sensing and Digital Earth at the Chinese Academy of Sciences for in situ data collection and data logging.
References
- 1. Wang MH. The Rayleigh lookup tables for the SeaWiFS data processing: accounting for the effects of ocean surface roughness. Int. J. Remote Sens. 2002;23(13): 2693–2702.
- 2. Tian LQ, Lu JZ, Chen XL, Yu ZF, Xiao JJ, Qiu F, et al. Atmospheric correction of HJ-1A/B CCD images over Chinese coastal waters using MODIS-Terra aerosol data. Sci. China Technol. Sci. 2010;53: 191–195.
- 3. Vanhellemont Q, Ruddick K. Advantages of high quality SWIR bands for ocean colour processing: Examples from Landsat-8. Remote Sens. Environ. 2015;161: 89–106.
- 4. Harrington JA, Schiebe FR, Nix JF. Remote sensing of Lake Chicot, Arkansas: Monitoring suspended sediments, turbidity, and Secchi depth with Landsat MSS data. Remote Sens. Environ. 1992;39: 15–27.
- 5. Chen J, Quan WT, Wen ZH, Cui TW. A simple ‘clear water’ atmospheric correction algorithm for Landsat-5 sensors. I: A spectral slope-based method. Int. J. Remote Sens. 2013;34(11): 3787–3802.
- 6. Nazeer M, Nichol JE, Yung YK. Evaluation of atmospheric correction models and Landsat surface reflectance product in an urban coastal environment. Int. J. Remote Sens. 2014;35(16): 6271–6291.
- 7. Mayo M, Gitelson A, Yacobi YZ, Ben-Avraham Z. Chlorophyll distribution in Lake Kinneret determined from Landsat Thematic Mapper data. Int. J. Remote Sens. 1995;16(1): 175–182.
- 8. Kloiber SM, Brezonik PL, Olmanson LG, Bauer ME. A procedure for regional lake water clarity assessment using Landsat multispectral data. Remote Sens. Environ. 2002;82: 38–47.
- 9. Wang DY, Feng XZ, Ma RH, Kang GD. A method for retrieving water-leaving radiance from Landsat TM image in Taihu Lake, East China. Chin. Geograph. Sci. 2007;17(4): 364–369.
- 10. Oyama Y, Matsushita B, Fukushima T, Matsushige K, Imai A. Application of spectral decomposition algorithm for mapping water quality in a turbid lake (Lake Kasumigaura, Japan) from Landsat TM data. ISPRS J. Photogramm. Remote Sens. 2009;64: 73–85.
- 11. Gong SQ, Huang JZ, Li YM, Wang HJ. Comparison of atmospheric correction algorithms for TM image in inland waters. Int. J. Remote Sens. 2008;29(8): 2199–2210.
- 12. Chen J, Fu J, Zhang MW. An atmospheric correction algorithm for Landsat/TM imagery basing on inverse distance spatial interpolation algorithm: A case study in Taihu Lake. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing. 2011;4(4): 882–889.
- 13.
Concha JA and Schott JR. In-water component retrieval over Case 2 water using Landsat 8: Initial results. 2014 IEEE Geoscience and Remote Sensing Symposium. 2014: 4458–4461.
- 14. Alcântara E., Bernardo N., Watanabe F., Rodrigues T., Rotta L., Carmo A., et al. Estimating the CDOM absorption coefficient in tropical inland waters using OLI/Landsat-8 images. Remote Sensing Letters. 2016; 7(7): 661–670.
- 15. Chen J, Zhu WN, Tian YQ, Yu Q. Estimation of colored dissolved organic matter from Landsat-8 imagery for complex inland water: case study of Lake Huron. IEEE Transactions on Geoscience and Remote Sensing, 2017; 55(4): 2201–2212.
- 16. Gordon HR, Wang M. Retrieval of water-leaving radiance and aerosol optical thickness over the oceans with SeaWiFS: A preliminary algorithm. Appl. Opt. 1994;33(3): 443–452. pmid:20862036
- 17. Li JS, Shen Q, Zhang B, et al. Retrieving total suspended matter in Lake Taihu from HJ-CCD near-infrared band data. Aquatic Ecosystem Health & Management. 2014;17(3): 280–289.
- 18. Wu CQ, Wang Q, Zhang YJ, et al. Remote sensing application system for water environments developed for Environment Satellite 1. Science China Earth Sciences. 2010;53(1): 45–50.
- 19.
Salama S, Monbaliu J. Atmospheric correction algorithm for CHRIS images application to CASI. International Symposium on Optical Science and Technology. International Society for Optics and Photonics. 2002).
- 20. Franz BA, Bailey SW, Kuring N, Werdell PJ. Ocean color measurements with the Operational Land Imager on Landsat-8: implementation and evaluation in SeaDAS. Journal of Applied Remote Sensing, 2015;9(1): 096070-1-096070-16.
- 21.
Vanhellemont Q, Ruddick K. Landsat-8 as a precursor to Sentinel-2: Observations of human impacts in coastal waters. ESA Special Publication, 2014: SP-726.
- 22.
Ruddick K. [Internet]. Ocean colour remote sensing in turbid coastal waters: CHL and TSM retrieval [cited 2017 Aug 27]. Available from: http://www.ioccg.org/training/SLS-2014/Ruddick_Lecture1_CHL&TSM_FINAL.pdf
- 23. Tang JW, Tian GL, Wany XY, Wang WM, Song QJ. The methods of water spectra measurement and analysis I: Above water method. J. Remote Sens. 2004;8(1): 37–44.
- 24. Mobley CD. Polarized reflectance and transmittance properties of wind-blown sea surfaces. Appl. Optics. 2015;54(15): 4828–4849.
- 25. Mobley CD. Estimation of the remote-sensing reflectance from above-surface measurements. Appl. Optics. 1999;38(36): 7442–7455.
- 26.
Solar Light Company, Inc. User's guide: Microtops II ozone monitor and sunphotometer (version 2.43). Doc. MTP05, Sol. Light Co., Glenside, PA, Oct. 2001.
- 27. Ångström A. The parameters of atmospheric turbidity. Tellus, 1964;16(1): 64–75.
- 28. Irons JR, Dwyer JL, Barsi JA. The next Landsat satellite: The Landsat Data Continuity Mission. Remote Sens. Environ. 2012;122: 11–21.
- 29. Pahlevan N, Lee ZP, Wei JW, Schaaf CB, Schott JR, Berk A. On-orbit radiometric characterization of OLI (Landsat-8) for applications in aquatic remote sensing. Remote Sens. Environ. 2014; 154: 272–284.
- 30. Kotchenova SY, Vermote E, Matarrese R, Klemm FJ. Validation of a vector version of the 6S radiative transfer code for atmospheric correction of satellite data. Part I: Path radiance. Appl. Opt. 2006;45(26): 6762–6774. pmid:16926910
- 31.
Vermote E, Tanre D, Deuze JL, Herman M, Morcrette JJ, Kotchenova SY. Second Simulation of a Satellite Signal in the Solar Spectrum-Vector (6SV), 6S User Guide Version 1. 2006.
- 32.
Zhang B, Li JS, Wang Q, Shen Q. Hyperspectral remote sensing of inland waters. Beijing: Science Press; 2012. p. 131–132.
- 33.
Gordon HR, Voss KJ. MODIS Normalized Water-leaving Radiance, Algorithm Theoretical Basis Document (MOD 18), version 4. NASA Contract Number NAS503163. 1999.
- 34. Pahlevan N, Lee ZP, Wei JW, Schaaf CB, Schott JR, Berk A. On-orbit radiometric characterization of OLI (Landsat-8) for applications in aquatic remote sensing. Remote Sens. Environ. 2014;154: 272–284.
- 35. Lv CG, Tian QJ, Wang L, Huang Y, Geng J. Estimation of Visible Surface Reflectance for Retrieving Aerosol Optical Depth Using Landsat-8 OLI Data. Remote Sens. Info. 2015;30(1): 43–50.
- 36.
landsat.usgs.gov [Internet]. Landsat Surface Reflectance Level-2 Science Data Products [cited 2017 Aug 27]. Differences between Landsat 4–7 and Landsat 8 Surface Reflectance algorithms. Available from: https://landsat.usgs.gov/landsat-surface-reflectance-data-products
- 37. Chen J, Cui TW, Lin CS. An improved SWIR atmospheric correction model: A cross-calibration-based model. IEEE Trans. Geosci. Remote Sens. 2013;52: 3959–3967.
- 38. Wang MH, Shi W, Tang JW. Water property monitoring and assessment for China's inland Lake Taihu from MODIS-Aqua measurements. Remote Sens. Environ. 2011;115(3): 841–854.
- 39. Chen J, Quan WT, Wen ZH, Cui TW. The use of MODIS 250 m bands to improve the MODIS 1 km ocean color atmospheric correction algorithm in turbid water. Advances in Space Research. 2012;51(9): 1750–1760.
- 40. Zhang MW, Ma RH, Li JS, Zhang B, Duan HT. A Validation Study of an Improved SWIR Iterative Atmospheric Correction Algorithm for MODIS-Aqua Measurements in Lake Taihu, China. IEEE Transactions on Geoscience and Remote Sensing. 2014;52(8): 4686–4695.
- 41. Gitelson A, Schalles JF, Hladik CM. Remote chlorophyll-a retrieval in turbid, productive estuaries: Chesapeake Bay case study. Remote Sens. Environ. 2007;109(4): 464–472.
- 42. Zhang F, Li J, Shen Q, Zhang B, Wu CQ, Wu YF, et al. Algorithms and schemes for chlorophyll a estimation by remote sensing and optical classification for Turbid Lake Taihu, China. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing 2015;8(1): 350–364.
- 43. Zhang B, Li JS, Shen Q, Chen D. A bio-optical model based method of estimating total suspended matter of Lake Taihu from near-infrared remote sensing reflectance. Environmental Monitoring and Assessment, 2008;145(1–3): 339–347. pmid:18066675