Improvements of 6S Look-Up-Table Based Surface Reflectance Employing Minimum Curvature Surface Method

We propose a methodology employing an interpolation technique on the Second Simulation of a Satellite Signal (6S) look-up table (LUT) to improve surface reflectance retrieval using Himawari-8/Advanced Himawari Imager (AHI). A minimum curvature surface (MCS) technique was used to refine the 6S LUT, and the solar zenith angle (SZA) and viewing zenith angle (VZA) increments were narrowed by 0.5°. The interpolation processing time was relatively short, about 3172 s per channel, and the interpolated xa and xb were well represented by the changes in SZA and VZA. An evaluation of the interpolated xa and xb for six cases revealed a relative mean absolute error of less than 5% for all channels and cases; however, a slight difference was evident for higher values of SZA and VZA. To evaluate the surface reflectance, we compared the surface reflectance derived using 6S LUT with that calculated using 6S only. Application of the interpolated 6S LUT showed a lower relative root mean square error (RRMSE) of 0.65% to 9.29% for all channels, than before interpolation. The improvement in surface reflectance measurements increased with the SZA. For a SZA above 75°, the RRMSE improved significantly for all channels (by 11.33–45.1%). In addition, when the MCS method was applied, the surface reflectance measurements improved without spatial discontinuity and showed good agreement with 6S results in a linear profile analyses. Thus, the method proposed can improve LUT based surface reflectance measurements in less time and increase the availability of surface reflectance data based on geostationary satellites.


Introduction
The spectral radiance measured using the satellite sensor is disturbed by absorption and scattering due to gases and aerosols present in the atmosphere on the path of sunlight from the Sun to the surfaces of the Earth and satellite-mounted sensors (Nunes et al. 2008;Hadjimitsis et al. 2010); this phenomenon is referred to as the "atmospheric effect". Absorption of sunlight by the atmosphere attenuates the radiation measured using the sensor. Scattering of sunlight can increase or decrease the radiation received by the sensor, depending on the atmospheric conditions (Proud et al. 2010). The magnitude of atmospheric effect is influenced not only by the amount of atmospheric constituents (i.e., the likelihood that a photon strikes an atmospheric component) but also by the atmospheric path length (i.e., the configuration of the Sun, Earth, and the satellite) (Hilker 2018). In satellite remote sensing, the atmospheric effect is a large source of uncertainty in observing the land surface. Accurate atmospheric correction is an essential pre-processing step to improve satellite image quality for land surface property analyses (Kaufman 1987;Jiménez-Muñoz et al. 2010).
To eliminate atmospheric effects, previous studies mainly applied the radiative transfer model (RTM), such as the Second Simulation of a Satellite Signal in the Solar Spectrum (6S) and MODerate resolution atmospheric TRANsmission (MODTRAN) (Adler-Golden et al. 1999; Responsible Editor: Soon-Il An. Kotchenova et al. 2006). The RTM method can be used to characterize the atmospheric effects of surface radiation signals measured using satellite sensors and has several advantages (Liang et al. 2010). It is not limited to a specific region or satellite sensor because it inputs the atmospheric conditions, geometry information, and sensor characteristics of areas where the atmospheric correction is to be performed. In addition, the RTM method is known to be accurate for simulating atmospheric effects in comparison with empirical line method (ELM) and improved dark-object subtraction (IDOS) (Zhou et al. 2011); thus, it is used in various fields, as well as in atmospheric correction of satellite imagery (Zhao et al. 2001;Liang et al. 2012;Zhang et al. 2016;Wu et al. 2017;Yeom et al. 2018). However, the RTM method has a disadvantage in that it takes a long time to perform atmospheric correction over a wide area, such as that of satellite images, because the calculation process is complex and the processing times are long. To overcome this issue, several studies have applied the look-uptable (LUT) method, which is an array-based method that replaces run-time calculations with simple array indexing (Liang et al. 2001;Lyapustin et al. 2011;Dorji and Fearns 2018). In constructing the LUT, the increment and range of input parameters directly affect the accuracy of atmospheric correction and computational efficiency (Qu et al. 2014).
The 6S LUT is used for the operational atmospheric correction algorithm of the Geostationary Operational Environmental Satellite-R series (GOES-R) and the Visual Infrared Imaging Radiometer Suite (VIIRS) (Liang et al. 2010;Franch et al. 2016). When constructing these LUTs, the increments for the solar zenith angle (SZA) and viewing zenith angle (VZA) were set to 4-5°, to optimize atmospheric correction accuracy and computational efficiency. At a higher SZA (>70°), the changes in atmospheric effects with variation in the angular measures (SZA and VZA) introduce significant error to surface reflectance measurements (Fraser et al. 1992). This is important not only in high-latitude regions where satellite observations are available at high SZA but also in terms of the availability of data measured with geostationary satellite sensors such as the Advanced Himawari Imager (AHI) on board the Himawari-8 satellite. With geostationary satellites in particular, observations are made at high SZA on a continuous basis; accurate surface reflectance data are required in terms of data utilization (Ruddick et al. 2014;Lee et al. 2015). To reduce the error at a high SZA, the interval of the input parameter should be set narrow in the LUT construction; however, this approach is limited by the long processing times.
In this study, we describe a method to refine the 6S LUT with short processing times that improve the quality of surface reflectance measurements over a wide SZA range, based on the minimum curvature surface (MCS) method. The MCS method is a common gridding technique used to construct a smooth surface from irregularly spaced data (Kaven et al. 2009;Rabah and Kaloop 2013). Computer implementation of the MCS method with a tension parameter to minimize extraneous inflection points was demonstrated by Smith and Wessel (1990); notably, the inclusion of a tension parameter did not introduce additional computational complexity to the MCS method. The MCS method is also commonly used for data interpolation in remote sensing fields (Patenaude et al. 2004;Nefeslioglu et al. 2012;Dietz et al. 2017;Gittings et al. 2017). Here, the MCS method was used to refine the 6S LUT to improve the surface reflectance quality for application to the Himawari-8/AHI that is a geostationary satellite observing the Asia-Oceania region over a wide range of SZAs and VZAs. Our results were compared to surface reflectance values calculated using 6S input values at the given observation time.
The remainder of this paper is organized as follows. Section 2 describes the data used in this research. The methodology for interpolating 6S LUT data is explained in Section 3. In Section 4, we present qualitative and quantitative assessments of our results, compared to surface reflectance calculations based on non-interpolated LUT and 6S. Our findings are summarized in Section 5.

Data
Himawari-8 is a next-generation geostationary meteorological satellite operated by the Japan Meteorological Agency (JMA) that was launched successfully on August 7, 2014 and began operation on July 7, 2015 (Yumimoto et al. 2016). The AHI on board the Himawari-8 satellite employs 16 channels from 0.47 to 13.3 μm (Choi and Suh 2018). The AHI instrument monitors a full disk area that includes Asia and Oceania; these areas are composed of various climatic zones and land cover types, ranging from broadleaf evergreen forests to barren lands. AHI images are acquired every 10 min, providing one of the highest temporal resolutions. Therefore, AHI provides multiple satellite images every day, which is a considerable advantage in areas where cloud cover is frequent, such as during Asia's summer season. AHI is suitable for the current study, because it acquires data on various atmospheric conditions, such as the water vapor content, ozone content, and aerosol optical depth (AOD), as well as measurements of the solar zenith and azimuth angles as the Earth rotates. For atmospheric correction, we used top of the atmosphere (TOA) radiance data from five shortwave channels centered at 0.47, 0.51, 0.67, 0.86, and 1.61 μm, respectively ( Table 1). The TOA radiance data used in this study has the same temporal resolution but different spatial resolution. The spatial resolution of two visible channels (channel 1 and 2) and one nearinfrared channel (channel 4) is 1 km. On the other hand, one visible channel (channel 3) and a SWIR channel (channel 5) have spatial resolutions of 0.5 km and 2 km, respectively. The TOA radiance data from the five shortwave channels were Korean Meteorological Society resampled at 2 km resolution. The spectral response function (SRF) for these channels, provided by the JMA, was also used for the 6S RTM simulations.
We performed atmospheric corrections of five shortwave channels only for cloud-free land surfaces. For the removal of cloud cover and water bodies, we used cloud mask and land/sea mask data provided by the National Metrological Satellite Center (NMSC). In addition, total precipitable water (TPW), total column ozone (TCO), and AOD data provided by the NMSC were used to characterize the atmospheric conditions at the AHI observation time; these data were generated by applying the GK-2A/AMI algorithm to Himawari-8/AHI. These data were used from 1 July to 7 August 2017 with 1 h intervals.
The cloud mask used in this research was output from the Geokompsat-2A (GK-2A)/Advanced Meteorological Image (AMI) cloud mask algorithm applied to Himawari-8/AHI. The GK-2A cloud mask algorithm was developed based on the threshold method considering channel characteristics according to the presence of the clouds and is currently used for operational purposes. Each pixel was classified by one of cloudy, probably cloudy and clear through threshold tests including single channel test, two-channel tests, spatial uniformity test and inversion layer correction (NMSC 2019). Its accuracy is 92% and 13% for POD and FAR, respectively, c om pa r e d t o t he M od e r a t e R e s o l u t i o n Im ag i n g Spectroradiometer (MODIS) cloud mask product (MYD35). The land/sea mask data was generated using MODIS yearly land cover data (MCD12Q1). Figure 1 shows the research flow chart for this study. We first constructed two LUTs (a non-interpolated and an interpolated 6S LUT) to remove atmospheric effects from AHI observation data. A non-interpolated 6S LUT was constructed for each AHI shortwave channel, through 6S RTM simulations of various input condition combinations. The interpolated 6S LUT was created by refining the interval of the input conditions in the non-interpolated 6S LUT, using the MCS method. The construction of the 6S LUT is discussed in more detail in the sections to follow. Then atmospheric correction coefficients of the most similar conditions to those of the present observation conditions were selected from two 6S LUTs. The surface reflectance values for the five shortwave channels were estimated using a simple equation provided by 6S. The surface reflectance is calculated only for the pixels classified as "confident clear" in the cloud mask product. To verify the improvement in the surface reflectance measurements, according to 6S LUT interpolation, the top-of-canopy surface reflectance calculated using the interpolated 6S LUT based on the MCS method (TOC reflectance MCS ) and the TOC surface reflectance calculated using the non-interpolated 6S LUT (TOC reflectance LUT ) were compared to the TOC surface reflectance calculated directly using 6S (TOC reflectance 6S ).

6S RTM
The 6S RTM was used to correct for atmospheric effects in TOA radiance measurements observed by AHI. The 6S was developed by Vermote et al. (1997) as an improved version of the 5S. A vector version of 6S (6SV) was introduced in 2005 and is based on successive orders of scattering approximations. In addition, 6SV is capable of accounting for radiation polarization based on calculations of the four Stokes vector components . The 6SV RTM calculates the scattering and absorption effects of the atmospheric components such as water vapor, ozone, and aerosol for a broad range of atmospheric, aerosol, spectral, and geometrical conditions. This calculation is performed by dividing the bandwidth of the channel by 2.5 nm. Kotchenova et al. (2008) evaluated the accuracy of four RTMs, including 6SV, MODTRAN, RT3 (radiative transfer), and SHARM (spherical harmonics), via a comparison with Coulson's tabulated values and Bréon's Monte Carlo code, both well-known benchmarks. 6SV shows higher simulation accuracy than MODTRAN, RT3 and SHARM with less than 1% error under several atmospheric condition scenarios. Due to the high accuracy of 6SV, it is used in radiometric calibration, retrieval of AOD, and surface albedo studies, as well as for atmospheric correction (Zhang et al. 2016;Wu et al. 2017;Yeom et al. 2018).
6SV calculates the surface reflectance based on the following: with where ρ TOA is the TOA reflectance. θ s , θ v and ϕ represent the SZA, VZA, and relative azimuth angle (RAA), respectively. T g represents the total gaseous transmission of H 2 O, CO 2 , O 2 and O 3 . S denote the spherical albedo of the atmosphere. ρ R + A represents the total reflectance of the molecule and aerosol scattering. T ↓ (θ s ) and T ↑ (θ v ) represent the value of the total transmission of the atmosphere along the radiation path between the Sun and the land surface and the satellite and land surface, respectively. Finally, ρ s is equation of transfer for a Lambertian homogeneous target of reflectance. 6SV provides three correction coefficients (xa, xb, and xc) as output for the input condition to remove atmospheric effects. The surface reflectance can be calculated from TOA radiance from Eq.
(3) using the three atmospheric correction coefficients.
In the above equation, ρ TOC denotes the calculated surface reflectance and L represents the TOA radiance (units: W m -2 μm -1 sr −1 ). xa, xb, and xc reflect the inverse of the transmittance, the scattering term of the atmosphere, and spherical albedo (atmospheric reflectance for isotropic light), respectively.

Pre-Construction of the 6S LUT
There are two methods for performing atmospheric correction using 6S. The first method is to perform atmospheric correction by calculating the atmospheric effects using 6S directly; the other approach is to simulate 6SV for various possible input conditions in advance and to store the results in the form of an LUT. The method of directly performing atmospheric correction using 6S is not suitable for processing a large amount of data, as similar calculations are repeated and the processing time is long. Therefore, we used the LUT method to avoid duplicate calculation processes and to accelerate the processing time. Table 2 shows the range and interval of input parameters used in this study for pre-construction of the 6S LUT; notably these values are important factors in atmospheric correction (Proud et al. 2010). The SZA and VZA had a range from 0°to 80°and varied in 5°increments. The RAA had a range of 0°to 180°, with a 10°i ncrement. For the atmospheric conditions, the ranges of TPW and TCO were set to 0-5 g cm −2 and 0.25-0.35 atm-cm, with intervals of 1 g cm −2 and 0.05 atm-cm, respectively, with reference to Lee et al. (2015). The input conditions for AOD consisted of 12 values an irregular interval, with reference to the values used by Liang et al. (2010). The three aerosol types were used: 'continental,' 'background desert,' and 'maritime'. We used the SRF of the AHI provided by the JMA to consider the spectral condition. The 6SV was performed for all combinations of each condition shown in Table 2. The three atmospheric correction coefficients (xa, xb, and xc) provided by 6S were stored in a LUT with a seven-dimensional structure answering input conditions (SZA, VZA, RAA, TPW, TCO, AOD, and aerosol type). To narrow the interval, we applied an MCS fitting method to the 6S LUT and narrowed the intervals for SZA and VZA conditions, as the volume of the constructed LUT becomes too large when the interval is refined for all input conditions. Among the three correction coefficients, xc did not change with SZA and VZA. Therefore, interpolation using MCS was not performed for xc. For the remaining correction factors xa and xb, we applied the MCS fitting method using an Interactive Data Language (IDL) routine to the 6S LUT to narrow the increment of SZA and VZA from 5°to 0.5°over the 80°range, taking into account the processing time and capacity. The MCS method is usually applied to define gradual variation in the surface and is equally applicable to equally spaced, line-oriented data, as a fast interpolation approach (Rufo et al. 2018). This method fits the linear elastic plate using data values that minimize the amount of bending to produce an optimal smooth surface (Grgić et al. 2015).

Accuracy Assessment
To evaluate the accuracy of surface reflectance according to the interpolation of the 6S LUT, we used three quantitative indicators, namely, mean absolute error (MAE), root mean square error (RMSE), and relative RMSE (RRMSE). The MAE is frequently used to evaluate accuracy and is derived from the absolute value of the difference between the LUT: look-up table, SZA: solar zenith angle, VZA: viewing zenith angle, RAA: relative azimuth angle, TPW: total precipitable water, TCO: total column ozone, AOD: aerosol optical depth, SRF: spectral response function The RMSE is a standard statistical metric for measuring accuracy (Chai and Draxler 2014) and is used to compare the differences between two values in various studies. The RRMSE is calculated by normalizing the RMSE using the mean value of the reference data and is expressed as a percentage (Despotovic et al. 2016). These indicators can be calculated by following equations: In the above equations, x represents the true value, in this case, the surface reflectance calculated directly from 6S; y is the estimated value, here, the surface reflectance derived by non-interpolated and interpolated 6S LUTs. N and x represent the number of observations and the mean of x, respectively. Figure 2 shows the xa and xb values of the interpolated 6S LUT and the non-interpolated 6S LUT, depending on SZA and VZA, for one combination of input parameters (relative azimuth angle, RAA: 30°; AOD: 0.1; TPW: 3 g/cm 2 ; TCO: 0.3 atm-cm, and aerosol type: continental). The non-interpolated 6S LUT constructed in this study well expressed the changes in the atmospheric correction coefficients with changes in the SZA and VZA. In the case of xa, the change in its value was more strongly affected by changes in the SZA than the VZA. In the case of xb, the effects of changes in the SZA and VZA on the change in xb were similar. When the SZA and VZA are relatively low, the rates of change in xa and xb are small in response to changes in the SZA and VZA. However, at higher SZA and VZA (>70°), the atmospheric correction coefficients change significantly with changes in the zenith angles. As mentioned, this tendency is an error factor in the atmospheric correction using the LUT. The interpolated 6S LUT shows a trend similar to that of the noninterpolated 6S LUT with the same minimum and maximum values. Particularly, xa and xb of the 6S LUT interpolated by the MCS fitting method well expressed the peak value at high SZA and VZA. We measured the elapsed time for interpolation of 6S LUT by the MCS method using a personal computer (Intel Core i5-9400F CPU; 16 GB memory). After constructing the 6S LUT, the time taken to construct the interpolated LUT by the MCS method was about 3172.18 s per channel. This is a much shorter time than that required when using 6S directly to construct an LUT with narrow intervals of SZA and VZA (0.5°). For quantitative evaluation of the interpolated 6S LUT, the atmospheric correction coefficients (xa and xb) from the interpolated 6S LUT were compared to the atmospheric correction coefficients calculated directly by 6S, as a reference for the six cases of the combination of the three RAA conditions (0, 90, 180°) and the two AOD conditions (0.01, 0.8). Table 3 shows the RAA and AOD conditions for each case. TPW and TCO were fixed at 1 g cm −2 and 0.3 atm-cm, respectively; these parameters have relatively little influence over atmospheric correction in the shortwave region, except for specific wavelengths, such as 1380 nm, that are absorbed strongly by water vapor (Meyer et al. 2016). Figure 3 shows the relative MAEs of xa and xb for each channel and case. The relative MAE is the value of the MAE divided by the average value of the reference data. In the case of xa, the interpolated xa through the MCS method showed good agreement with the xa calculated directly using 6S in all cases and for all channels, with a relative MAE of less than 2%. The minimum value of the relative MAE was 0.478% in case 1 of channel 4, and the maximum value was 1.4% in case 4 of channel 5. Although the differences were small under most SZA and VZA conditions, the difference in xa became increasingly large when the SZA exceeded 75°, up to 0.014 (8.15%) in case 6 of channel 5. For SZA above 75°, the xa increased sharply with the SZA. The interpolation performance of the MCS method becomes degraded when the data change over short intervals (Rufo  ; thus, xa of the interpolated 6S LUT differed slightly from xa calculated directly by 6S. As the condition of RAA changes in all channels, the relative MAE value of xa fluctuates within 0.01%. However, as the AOD condition changes, there is a relatively large error variation (0.5%). The reason for this error tendency is that xa is more sensitive to AOD changes than RAA changes. xb also showed a relative MAE of less than 5% in all cases and channels, compared to xb calculated directly with 6S. In all cases, a large relative MAE was observed in channel 5, up to 2.8% for case 6. Given the same AOD condition for all channels, a small relative MAE was obtained for RAA equal to 90°(cases 3 and 4). The relative MAE was 0.18% in case 4 of channel 1. In the case of xb, similar to xa, the differences were small for most SZA and VZA conditions; however, the difference increased significantly when SZA and VZA exceeded 75°, up to −1.143 (−16.8%) in case 6 of channel 5. In the case of SZA and VZA above 75°, xb increased sharply with SZA and VZA. Therefore, due to the limitations of MCS interpolation mentioned above, the xb of the interpolated 6S LUT differed slightly from xb calculated directly by 6S under high VZA and SZA conditions. The relative MAE of xb has variability as the RAA and AOD change. Thus, xb is more sensitive to changes in RAA than xa.

Evaluation of Accuracy for Surface Reflectance
Calculating surface reflectance directly with 6S can take a long time. Therefore, we quantitatively evaluated the LUT based surface reflectance by sampling pixels corresponding to various SZA and VZA conditions. We sampled about 58,000 pixels per channel for five shortwave channels of the AHI. Figure 4 shows the VZA and SZA distributions of sampled pixels for quantitative evaluation of the surface reflectance. The sampled pixels were relatively uniform in number for all SZAs, with the exception of SZA < 15°. For the VZA, however, the number of samples was small for VZA < 20°and VZA > 55°. AHI is located at 140.7°above the equator for its Earth's observations. Therefore, in the region corresponding to a low VZA range, there is a small number of land pixels; in the region corresponding to a high VZA range, the land pixel is relatively small. In particular, equatorial regions with low VZA show low sampling numbers due to the large number of pixels contaminated by clouds. Table 4 shows MAE, RMSE, and RRMSE for TOC reflectance LUT and TOC reflectance MCS compared to TOC reflectance6S for all sampled data. The MAE and RMSE of TOC reflectance LUT range from 0.0031-0.0050 and from 0.0043-0.0065, respectively. The MAE and RMSE of TOC reflectance MCS were 0.0014-0.0026 and 0.00194-0.0036, respectively, showing higher accuracy after 6S LUT interpolation of all five channels. The RRMSE of the 6S LUT varied from 2.10% to 17.51% before interpolation and from 0.65% to 9.29% after interpolation. After interpolation, the RRMSE of all channels improved by approximately two-fold, compared to that before interpolation. The greatest improvement in RRMSE was 7.92% on channel 1. Channels 4 and 5 showed very small RRMSE values of less than 1%.
Because surface reflectance is mainly affected by the SZA and VZA, we evaluated the accuracy of surface reflectance calculated using 6S LUT according to the SZA conditions. For this evaluation, the SZA was categorized into 17 classes in 5°intervals. Figure 5 shows the MAE, RMSE, and RRMSE of TOC reflectance LUT and TOC reflectance MCS for five channels, according to the 17 SZA classes. In the case of TOC reflectance LUT , the MAE, RMSE, and RRMSE tended to increase with the SZA for all channels; this increase was particularly sharp for high SZA values, in which SZA > 60°. The RMSE values of channels 1-5 ranged from 0.00264-0.0606, 0.00206-0.0466, 0.000138-0.0336, 0.00172-0.0469, and 0.00134-0.0356, respectively. RRMSE also varied greatly according to the SZA, and increased by 72.6%, 64.9%, 30.4%, 16.2%, and 15.1% in the largest SZA class for channels 1-5, respectively. These findings show that the difference between TOC reflectance LUT and TOC reflectance 6S increased with the SZA and showed large differences under high SZA conditions. For TOC reflectance MCS , the MAE, RMSE, and RRMSE were lower than those of the TOC reflectance LUT for all SZA classes and all channels. The MAE, RMSE, and RRMSE also increased with the SZA; however, with the exception of the largest SZA class, there was no significant increase in the error with SZA, unlike that observed for the TOC reflectance LUT . The RMSE values for the largest SZA class were 0.03, 0.0182, 0.0171, 0.0105, and 0.0088, respectively, for channels 1-5, which are relatively larger than other SZA classes. As mentioned above, when the SZA varies between 75°and 80°, the atmospheric correction coefficient changes significantly according to the SZA, and a relatively large error occurs when interpolating the 6S LUT with the MCS method. Despite these errors, the TOC reflectance MCS showed a significantly lower RRMSE than the RRMSE of TOCreflectance LUT by 45.1%, 39.69%, 19.8%, 12.6%, and 11.3%, respectively, in the highest SZA class.
The improvement in the accuracy of the surface reflectance was due to the improvement in the xa and xb values, using the 6S LUT interpolation. Figure 6 shows the accuracy indices for xa and xb of channels 1 and 2, according to the SZA class. The accuracy of xa displayed a trend similar to that of TOC reflectance. As the SZA increased, the accuracy decreased; however, the accuracy was much higher after 6S LUT interpolation. In channels 1 and 2, the RRMSE of xa had a low error range of 1.30% (first SZA class) to 3.99% (last SZA class) and 1.24% (second SZA class) to 4.23% (last SZA class), respectively, (a) (b) Fig. 7 Spatial distribution of (a) TOC reflectance LUT and (b) TOC reflectance MCS at channel 3 on August 4, 2017 at 07:00 UTC. The grey solid contour line corresponds to SZA in degree units. TOC reflectanceLUT: top-of-canopy surface reflectance calculated using the non-interpolated 6S LUT; TOC reflectanceMCS: TOC surface reflectance calculated using the interpolated 6S LUT, based on the minimum curvature surface method Korean Meteorological Society depending on the SZA class. After interpolating the 6S LUT, the accuracy of xb of all SZA classes improved, and the RRMSE decreased by 13.3% and 14.1% for channels 1 and 2 respectively, in the last SZA class. However, the improvement in xb was less than that of xa. Because xa is more sensitive to the SZA and VZA than other input parameters, it is greatly improved by interpolation of the 6S LUT for SZA and VZA. On the other hand, given that xb is sensitive to both SZA and VZA, as well as RAA and AOD, the difference in xb can be relatively large, depending on the difference between the LUT construction condition and the actual condition. For this reason, relatively low improvement was achieved in xb than in xa. Therefore, the improvement in TOC reflectance by interpolation of the 6S LUT was mainly due to the improvement in xa. Similar trends were observed in channels 3-5; however, they were not included for brevity. Fig. 8 Linear profile of surface reflectance for channels (a) 3, (b) 4, and (c) 5. The red, blue, green circle at the top of each figure indicates TOC reflectance LUT , TOC reflectance MCS , and TOC reflectance 6S of each channel, respectively; red and blue triangles at the bottom of each figure indicate the difference between TOC reflectance LUT and TOC reflectance 6S , and the difference between TOC reflectance MCS and TOC reflectance 6S , respectively. TOC reflectance6S: TOC reflectance calculated directly using 6S Korean Meteorological Society The spatial distribution of the surface reflectance calculated from the Himwari-8/AHI TOA data was analyzed to confirm the improvement with interpolation of the 6S LUT. The TOC reflectance LUT and TOC reflectance MCS for channel 3 on August 4, 2017 at 07:00 UTC are shown in Fig. 7. For the TOC reflectance LUT , a discontinuity occurred in the spatial distribution of the surface reflectance in Australia. This is similar to the SZA distribution. The larger the SZA, the larger the difference between neighboring pixels. The discontinuity was clearly visible at higher SZA above 65°, according to the interval of SZA (5°) used for construction of 6S LUT in this study. As mentioned above, the larger the SZA, the greater the change in the atmospheric correction coefficients as the SZA changes. Therefore, this discontinuity appears to be caused by the difference between the LUT construction condition and the sensor observation condition. On the other hand, the TOC reflectance MCS did not show discontinuity in the surface reflectance with SZA variation, particularly for a high SZA. For this reason, there is a difference between the TOC reflectance MCS and TOC reflectance LUT . Other channels also show the same spatial distribution but are excluded for brevity. Figure 8 shows the linear profiles of the TOC reflectance LUT , TOC reflectance MCS , and TOC reflectance 6S for channels 3-5 for the region marked with an orange solid line in Fig. 7, where discontinuity occurs in the spatial distribution of surface reflectance. In Fig. 8, the red, blue, and green circles correspond to TOC reflectance LUT , TOC reflectance MCS , and TOC reflectance 6S , respectively; red and blue triangles represent the differences between TOC reflectance LUT and TOC reflectance 6S and between TOC reflectance MCS and TOC reflectance 6S , respectively. TOC reflectance LUT and TOC reflectance MCS showed a linear profile similar to TOC reflectance 6S as a whole; however, there was a difference in TOC reflectance LUT as the SZA at the time of observation differed considerably from that of the SZA used for LUT construction. When the actual SZA is smaller (larger) than the SZA corresponding to the LUT condition, the difference indicates a positive (negative) value. The magnitude of the difference increased with the SZA; the RRMSE values of channels 3-5 were 14.16%, 11.03%, and 10.06%, respectively, for a SZA above 70°. TOC reflectance MCS showed a linear profile similar to the TOC reflectance 6S for SZA ≤ 75°. For a SZA above 75°, TOC reflectance MCS tended to be overestimated compared to the TOC reflectance 6S ; however, each channel showed a low MAE of 0.007, 0.006, and 0.007, respectively. In the linear profile of TOC reflectance MCS , the difference between TOC reflectance MCS and TOC reflectance 6S also repeated the positive/negative trend depending on whether the SZA value (interval chosen) corresponded to a small or large SZA; this is an inherent error of the LUT method. However, the difference between TOC reflectance MCS and TOC reflectance 6S was small.
For the TOC reflectance MCS with SZA > 70°, the RRMSEs of channels 3, 4, and 5 were 3.55%, 2.01%, and 1.84%, respectively, which improved by 10.16%, 9.02%, and 8.22%, respectively, compared to the RRMSE of the same channel before interpolation of the 6S LUT. Similar trends were observed in channels 1 and 2; however, they were not included for brevity.

Concluding Remarks
We proposed a method for improving the quality of surface reflectance by interpolating the 6S LUT for application to the Himawari-8/ AHI. The interpolated 6S LUT showed lower MAE in all channels compared to the non-interpolated 6S LUT. TOC reflectance MCS more agreed well with TOC reflectance 6S in all cases than TOC reflectance LUT . Particularly, for SZA > 75°, an improvement of 45.1%, 39.6%, 19.8%, 12.57%, and 11.33% in RRMSE were indicated for channels 1-5, respectively. Linear profile analyses indicated a slight overestimation of TOC reflectance MCS for SZA > 75°, but showed good agreement without discontinuity along 6S LUT interval over the entire SZA range.
The proposed method in this study can improve the accuracy of surface reflectance measurements, even for a high SZA, by constructing a detailed 6S LUT with a short processing time. This can increase the availability of surface reflectance values at high latitudes, where high SZA observations are made. In addition, the proposed method can be applied to next-generation geostationary orbiting satellites, such as GK-2A (South Korea), GOES-R (United States), and Meteosat Third Generation (Europe). In order to apply this method to next-generation geostationary satellites, 6S LUT for each satellite must be built in advance, which takes a long time. However, if the 6S LUT is constructed, the proposed method can be applied without any optimization method and the production of improved surface reflectance can be expected.
In order to improve the accuracy of the calculation of the surface reflectance in future studies, the accuracy of the atmospheric correction coefficient should be ensured for the portion where the value of the atmospheric correction coefficient varies greatly with SZA and VZA. In addition, it is necessary to analyze and improve the error according to the LUT interval of other input variables such as AOD and RAA, which greatly affect xb.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.