Are visibility-derived AOT estimates suitable for parameterizing satellite data atmospheric correction algorithms?

Aerosol optical thickness (AOT) is an important parameter in radiative transfer models (RTMs) used for atmospheric correction of remotely-sensed data. It is often estimated from horizontal visibility measurements by use of the Koschmieder formula or other related methods built into RTMs. This article assesses the accuracy of this estimation, in the context of atmospheric correction, by comparing AERONET AOT data with AOT estimates from UK Met Office visibility data at a site in Hampshire, UK. Root mean square errors are calculated for a number of visibility categories (0–10, 10–20, 20–30, and 30–40 km) and are found to be high for all visibilities (ranging from half to more than double the mean AOT for each category). For all visibilities km, these errors are significantly higher than those from other AOT estimation methods. Simulations performed with the 6S RTM show that the effect of these errors on satellite-level radiances are large (up to 36 ), and the change in vegetation indices (normalized difference vegetation index and atmospherically resistant vegetation index) is smaller, but still significant. It is recommended that estimations of AOT based upon visibility measurements are only used if no alternatives are available, and that great caution is used when estimation is performed for visibilities km.


Introduction
Atmospheric correction is an essential part of many remote-sensing studies, as it removes the perturbing effect of the atmosphere from the data and allows the measurements to be determined accurately in physical units (Slater 1980). This is essential as remotely-sensed data are now used to generate many quantitative products, such as those in the Global Climate Observing System Essential Climate Variables (GCOS 2004). Many atmospheric correction methods are available, including image-based approaches such as dark object subtraction (Chavez 1988), empirical relationships between sensor and ground reflectances such as the empirical line method (Smith and Milton 1999), and simulations based on radiative transfer models (RTMs) such as 6S (Vermote et al. 1997) or Moderate Resolution Atmospheric Transmission (MODTRAN) (Berk et al. 1999). These physical methods use a RTM to simulate the effect of the atmosphere on the light passing through it, and require parameterizing to simulate the atmospheric conditions at the time of image acquisition, with the accuracy of this parameterization controlling the accuracy of the resulting correction.
A key parameter is the aerosol optical thickness (AOT) (a measure of the reduction in atmospheric clarity due to the presence of aerosols), but AOT measurements coincident with image acquisition are often unavailable, so AOT is instead estimated from horizontal visibility measurements. There is uncertainty about the relationship between horizontal visibility and AOT, and this article investigates the accuracy of these AOT estimates and the atmospheric correction errors which may result from using an incorrect AOT parameterization. This is particularly important now that more data are being automatically atmospherically corrected before release by satellite operators.

Background
2.1. Why are visibility measurements used? AOT is measured accurately, and automatically, by a global network of over 300 Cimel CE-318 sun photometers (Holben et al. 1998). These accurate measurements are used for RTM parameterization where possible, but many satellite images do not include an AERONET site location. Only 6.7% of Landsat scene footprints over land (8.5% excluding Antarctica) have ever had an AERONET site located within them, and furthermore, the period of record at many AERONET sites is short, thus making it difficult to parameterize atmospheric corrections for archive imagery, which stretches back to the early 1970s for some satellites. AOT measurements can be taken from handheld instruments such as the Microtops II sun photometer (Morys et al. 2001); however, this requires researchers to be in the field at the time of satellite image acquisition which is often challenging and impossible for archive imagery.
In comparison, measurements of horizontal visibility are acquired automatically at the majority of airports worldwide and at automated meteorological stations run by national and international meteorological agencies. Measurements are often acquired very frequentlyup to every minute, compared with every 15 min for AERONET dataand there are long periods of record at many sites. In general, it is far more likely that, for a given satellite image, coincident visibility data will be available than AERONET-based AOT data ( Figure 1).

How is visibility measured?
The word 'visibility' has a wide range of meanings, both in everyday speech and in scientific research. When measured and recorded for meteorological purposes, one of many specific technical definitions of visibility applies. The observer visibility, V 0 , is the greatest distance at which it is just possible to see and identify, with the unaided eye, a prominent dark object against the sky at the horizon. An alternative, more specifically defined measure, is the visual range which is defined as 'the distance, under daylight conditions, at which the apparent contrast between a specified type of target and its background becomes just equal to the threshold contrast of an observer' (American Meteorological Society 2013). This is dependent on the minimum contrast that the observer can discriminate at the moment of observation, which varies significantly between people, and thus reduces the comparability of these measurements.
The meteorological range, V, is defined using the following formula (Koschmieder 1925), where σ is the horizontal extinction coefficient, is the threshold contrast, and C 0 is the inherent contrast of the object. Koschmieder suggested that under most viewing conditions, the threshold contrast is 0.02, and proposed a theoretical 'black object' whose inherent contrast is −1. Substituting these gives the standard Koschmieder formula, Thus, V is directly related to σ and is usually recorded by automated visiometers which directly measure σ. Use of the meteorological range removes the subjectivity of human observation, but the resulting visibility is dependent on the parameter assumptions made by Koschmieder. Although Koschmieder suggested a threshold contrast of 0.02, there have been suggestions that the appropriate value could range from 0.01 to 0.2 (World Meteorological Organization 2010).
The majority of modern visibility measurements are recorded as the meteorological optical range (MOR), which is defined as 'the path length in the atmosphere required to reduce the luminous flux in a collimated beam from an incandescent lamp, at a colour temperature of 2700 K, to 5 per cent of its original value' (World Meteorological Organization 2010, 1.9-2). This is equivalent to the meteorological range with a threshold contrast of 0.05 (World Meteorological Organization 2010), In general, observer visibility, V 0 , is lower than the meteorological range, V, with Carr (2005) developing an approximate relationship between them: As MOR is the standard meteorological measure of visibility, all mentions of visibility in the rest of this article should be taken to mean MOR.

How is AOT estimated from visibility?
AOT is defined as the integrated vertical extinction coefficient over a vertical atmospheric column of unit cross section. Thus, to accurately derive AOT, the vertical profile of the extinction coefficient must be known, and a number of relationships between AOT and visibility which take into account various measures of the vertical extinction profile, such as the planetary boundary layer height, have been derived in recent years (Kessner et al. 2013). However, none of these methods have been implemented in atmospheric correction software, as thus the rest of this article will focus on the simpler relationships which tend to be used in atmospheric correction.
The simplest relationship between AOT and visibility is also based on the Koschmieder equation (Equation 2), with a replacement of σ, the horizontal extinction coefficient, with AOT (τ), This simple substitution of terms assumes that (1) the horizontal and vertical extinction coefficients are equal, (2) the relationship between horizontal visibility and horizontal extinction coefficient also holds for horizontal visibility and the vertical extinction coefficient, and (3) the vertical extinction profile has no effect on the AOT value. Chan (2009) found that these assumptions did not hold for data collected over Hong Kong International Airport during winter 2008-2009. A number of occasions were observed when the vertical extinction coefficient increased significantly, but the horizontal extinction coefficient changed very little. Measurements with a ground-based Doppler lidar showed that this was caused by a significant change in the vertical extinction profile, with high aerosol concentrations at altitudes of 2-4 km. These increased concentrations were not associated with any noticeable atmospheric phenomena such as clouds or dust storms, and thus could not be observed through visibility measurements. Similar situations can be caused by other unexpected layering in the atmosphere, such as temperature inversions.
Modified versions of the Koschmieder equation have been proposed, using a variety of threshold contrast values, or adding a correction factor. For example, Vermote (2002), cited in Retalis et al. (2010) added a correction factor ðþ0:08498Þ to the standard equation.
Two widely used RTMs, 6S (Vermote et al. 1997) and MODTRAN (Berk et al. 1999), calculate AOT by numerically integrating extinction coefficients for many levels of the atmosphere. 6S interpolates the particle density at 32 levels in the atmosphere from two profiles taken at 5 and 23 km visibility, before converting to extinction coefficient and then integrating to get AOT (Vermote et al. 2006). MODTRAN also calculates AOT as the integral of extinction coefficients, using data at 35 levels. These extinction coefficients are produced by scaling the values in the standard aerosol profiles. The Koschmieder formula with a correction for Rayleigh scattering, is used to calculate the scaling factor at ground level, S 0 , and the rest of the scaling factors are dependent on this and the season and volcanic conditions (Carr 2005).

Previous work
A number of previous studies have compared AOT measurement and visibility-based estimates of AOT, although relatively few studies have been performed with a focus on using visibility-derived AOT estimates for atmospherically corrected satellite images. Regardless of the application of the resulting AOT data, the justification for using visibility-based estimates was the ease of data availability and the long time-series of data available (1923 to present for a number of stations; NCDC 1998). Retalis et al. (2010) found a significant relationship between AOT from a Microtops sun photometer and AOT estimated from visibility using the Vermote et al. (1997) formula (r ¼ 0:76), but this relationship relied heavily on a few extreme AOT values. So, Cheng, and Tsui (2005) found a weaker correlation (r ¼ 0:38) between Moderate Resolution Imaging Spectroradiometer (MODIS) AOT data and visibility-based AOT estimates and suggested a number of reasons for the poor relationship, including the low spatial resolution of the MODIS data and the assumptions detailed earlier. A higher correlation (r ¼ 0:89) was found by Chen et al. (2009) between MODIS AOT data and estimated AOT based on the Koschmieder formula with corrections for atmospheric pressure and relative humidity. These previous studies examined the correlation between measured and estimated AOT values, whereas this study will quantify the error in real-world units between measurements collected at the same time and location, thus giving a better assessment of the accuracy of the estimation methods.

Data sources
Spatially coincident measurements of AOT and visibility are recorded at the STFC Chilbolton Facility for Atmospheric and Radio Research (CFARR) in Hampshire, UK (51.1457N, 1.4284W). Furthermore, a meteorological station run by the UK Met Office (UKMO) also recording visibility data is located approximately 10 km away at Middle Wallop (51.1498N, 1.5699W).

AOT measurements
AOT at a number of wavelengths was measured and logged every 15 min at CFARR as part of NASA's AERONET network (Holben et al. 1998). We include only cloud-free data (AERONET Level 2 data) in this study as the presence of clouds invalidates the assumptions of the visibility-AOT relationship. The accuracy of the AOT recorded by the Cimel instrument is AE 0:01 (Eck et al. 1999), and it is automatically corrected for Rayleigh scattering.

CFARR visibility measurements
Visibility data, measured as MOR, were logged every minute at CFARR from September 2009 until October 2011, with over 1 million individual measurements (Science and Technology Facilities Council 2011). The visibility instrument (Vaisala PWD21,Vaisala,Helsinki,Finland) can only measure a maximum visibility of 20 km; thus, on clear days, the instrument saturates, making comparisons with AERONET data difficult. There were only 45 days when the visibility sensor did not saturate during the period of data availability. The visibility measurements are accurate to approximately AE10% from 0 to 10 km and AE15% from 10 to 20 km (Vaisala 2010), within the 10-20% considered to be a realistically achievable accuracy (Crosby 2003).

6S and MODTRAN simulations
The original FORTRAN code for 6S is freely available, and the relevant functions for estimating AOT from the input visibility were rewritten in R and incorporated into the analysis code.
The code for MODTRAN is not available, and so simulations were performed at many visibilities (0-50 km in steps of 0.5 km) and the results linearly interpolated to provide AOT estimates for every measured visibility. As MODTRAN does not output the estimated AOT, it was calculated by numerically integrating the extinction coefficient at each altitude and subtracting the Rayleigh optical depth, τ R , which was calculated from (Dutton et al. 1994), τ R ¼ P 1013:25 0:00877λ À4:05 ; using the MODTRAN-provided ground pressure value, P, at a wavelength, λ, of 500 nm.

Validation method
Time series of each data source were analysed using the zoo (Zeileis and Grothendieck 2005) and xts (Ryan and Ulrich 2011) packages in the R statistical programming environment (R Development Core Team 2011) to select and merge all AOT and visibility measurements within 15 min of each other and then produce error statistics. All of the codes used to produce the results in this article are available at https://github.com/robintw/ VisAOT, and the ProjectTemplate framework (White 2012) has been used to improve reproducibility. Overall, 4677 matched measurements were found, across the period from October 2005 until July 2013.

Simulation of resulting atmospheric correction errors
Radiative transfer modelling was used to assess the effect of this error on Landsat atsensor radiances and vegetation index values. A simulated pixel with a typical vegetation reflectance spectrum was used, and simulations were performed using the Py6S interface (Wilson 2012) to the 6S model (Vermote et al. 1997) using the parameters shown in  (Kaufman and Tanré 1992) with the ARVI γ value set to 1.0, recommended as the optimum value by Kaufman and Tanré (1992).

Results and discussion
Comparing the three AOT estimation methods (Figure 2) shows the expected exponential trend for all methods. At low visibilities <10 km, the standard Koschmieder equation and 6S produce very similar results, but results from MODTRAN are significantly higher. The situation is different for high AOTs, with 6S producing the highest AOTs and MODTRAN the lowest.
The visibility measurements at CFARR and UKMO Middle Wallop are comparable (r ¼ 0:78 and very similar qualitative trends, Figure 3), albeit with a slight temporal mismatch which varies in magnitude and direction. Therefore, UKMO data were used as the wider range of visibilities (up to 40 km) includes the clear conditions during which the majority of remote-sensing measurements are likely to be collected.
Error between the estimated AOT from visibility data, calculated using each of the estimation methods, and the AERONET data was calculated ( Table 2). The m values show whether the method overestimates (m > 1) or underestimates (m < 1). The visibility data were categorized to allow the analysis of the error at different visibilities, and Figure 4 shows examples of webcam images indicating different visibilities, to allow these categorizations to be related to real-world conditions.
The results show an inverse relationship between visibility and RMSE, but even in the highest visibility category, the RMSE is still over 50% of the mean AOT for that category. The RMSEs calculated over the whole visibility range show that MODTRAN performs worse than Koschmieder and 6S, which produce a similar error. For all except the lowest visibility category, the Koschmieder formula has the lowest RMSE. All methods overestimate significantly for low visibilities, with RMSEs almost twice the average AOT value and with MODTRAN producing significantly worse results than 6S and the Koschmieder formula. This general poor performance is likely to be caused by failure of the underlying assumptions regarding the equality of horizontal and vertical extinction profiles.
In the range of visibilities where remote-sensing data acquisition is most likely to take place, >20 km, the RMSEs are still a large proportion of the mean AOT value and the Koschmieder formula and MODTRAN tend to underestimate, whereas 6S overestimates.
Other AOT estimation methods also have significant errors. For example, the LEDAPS method for obtaining AOT from Landsat images produces results in which approximately 56% of measurements are within AE0:1 of co-located AERONET measurements (Maiersperger et al. 2013). The MODIS MOD04 product also estimates AOT from the image data itself, although at a lower resolution, with 66% of values within 0:05 AE 15% (Chu et al. 2002), which corresponds to a representative error for a moderate AOT (0.2) of AE 0:08. AERONET measurements are the most accurate, with an error of approximately AE 0:01 (Eck et al. 1999). These are approximately equivalent to the errors in visibility-derived AOT when the visibility is >30 km, but for lower visibilities, the error in visibility-derived AOT is significantly higher than the errors from other methods.
Modelling was used to assess the effect of this error on at-sensor radiances and the resulting NDVI values. A simulated pixel was configured with the built-in standard green vegetation spectrum in the 6S model, with other parameters set as in Table 1. Simulations were performed using the average AOT for each of the visibility ranges, plus this AOT perturbed both positively and negatively by the average AOT error across all estimation methods for that visibility range. The resulting at-sensor radiances, in the Landsat ETM+ bands, and NDVI and ARVI values calculated using these radiances are shown in Table 3.
The results show that at-sensor radiances can be significantly affected by AOT perturbations, with larger errors at shorter wavelengths. Comparing Table 3 to the Noise Equivalent Delta Radiance (NEΔL) values shown in Table 4 shows that the radiance error is always greater than the NEΔL, with a maximum error of over 30 times the NEΔL, and an error at moderate visibilities of 3-4 times the NEΔL. As expected, the vegetation indices are affected less than the radiances, but there are still significant errors of up to 0.12 in NDVI and 0.09 in ARVI, with errors for moderate visibilities of approximately 0.02. These errors are relatively small, but even small differences may affect estimates of biomass production, for example, Kaufman and Holben (1993) found that a NDVI difference of 0.04 corresponded to biomass production errors of 11-30%. ARVI generally has lower errors than NDVI, but the difference is smaller than expected, and shows that ARVI is still susceptible to significant atmospheric effects.

Conclusions
This study suggests that the estimation of AOT from horizontal visibility measurements can produce significant errors. This is based upon the data set examined in this study, and so errors will vary in different situations. However, we believe that this data set is representative of the AOT-visibility relationship under standard mid-latitude atmospheric conditions, and thus these errors are likely to be representative of the magnitude of potential errors in many atmospheric correction applications. Overall, errors are highest at lowest visibilities, but even at high visibilities, there can be errors of 50% of the AOT value. In all but the lowest visibility category (0-10 km), the raw Koschmieder formula had the lowest error, followed by the 6S method, with the MODTRAN method producing the highest error (particularly at low visibilities where the Notes: The AOT error used is the mean of the AOT errors from each AOT-derived visibility estimation method (from the final column in Table 2). τ is the AOT used for the simulation (set to zero if the perturbation would result in a negative AOT), with Δτ showing the AOT perturbation used, are the radiances at satellite level for the Landsat ETM bands, and the NDVI and ARVI are also shown. error from MODTRAN is significantly higher than the other methods). Care should be taken when comparing results producing using different methods, as the differences between estimation methods mean that simulations carried out in 6S and MODTRAN using the same visibility may produce significantly different results due to the differing estimates of AOT. Remote-sensing data are not normally acquired during conditions of very low visibility (e.g. the first two images in Figure 4); therefore, the errors seen at 0-10 km visibilities are unlikely to be experienced with operational remote sensing. However, at the moderate visibilities in which much remote-sensing acquisition is performed, the error in the AOT is around 70%. For visibilities <30 km, the errors in visibility-derived AOT are significantly higher than those from alternative AOT estimation methods (including image-based methods and standard satellite products), and for higher visibilities, the errors are approximately equal. Thus, other methods for estimating AOT should be used where possible, with visibility-derived AOT only used as a last resort. Furthermore, great care should be taken when using visibility-derived AOT at low visibilities.
Simulated radiance values performed using the AOT errors show that the resulting radiance errors are significantly larger than the NEΔL, with errors of up to 36 Wm À2 sr À1 . Even at the highest visibilities, the radiance error is twice the NEΔL. Errors in vegetation indices are lower, but still significant for both NDVI and ARVI. Thus, if it is suspected that inaccurate AOTs have been used to atmospherically correct an image, radiance values and vegetation indices should be used with care.