Novel ray tracing method for stray light suppression from ocean remote sensing measurements

We developed a new integrated ray tracing (IRT) technique to analyze the stray light effect in remotely sensed images. Images acquired with the Geostationary Ocean Color Imager show a radiance level discrepancy at the slot boundary, which is suspected to be a stray light effect. To determine its cause, we developed and adjusted a novel in-orbit stray light analysis method, which consists of three simulated phases (source, target, and instrument). Each phase simulation was performed in a way that used ray information generated from the Sun and reaching the instrument detector plane efficiently. This simulation scheme enabled the construction of the real environment from the remote sensing data, with a focus on realistic phenomena. In the results, even in a cloud-free environment, a background stray light pattern was identified at the bottom of each slot. Variations in the stray light effect and its pattern according to bright target movement were simulated, with a maximum stray light ratio of 8.5841% in band 2 images. To verify the proposed method and simulation results, we compared the results with the real acquired remotely sensed image. In addition, after correcting for abnormal phenomena in specific cases, we confirmed that the stray light ratio decreased from 2.38% to 1.02% in a band 6 case, and from 1.09% to 0.35% in a band 8 case. IRTbased stray light analysis enabled clear determination of the stray light path and candidates in in-orbit circumstances, and the correction process aided recovery of the radiometric discrepancy. ©2016 Optical Society of America OCIS codes: (010.0280) Remote sensing and sensors; (290.2648) Stray light. References and links 1. Q. Guo, J. M. Xu, and W. J. Zhang, “Stray light modelling and analysis for the FY‐2 meteorological satellite,” Int. J. Remote Sens. 26(13), 2817–2830 (2005). 2. H. Tonooka, “Inflight straylight analysis for ASTER thermal infrared bands,” IEEE Trans. Geosci. Rem. Sens. 43(12), 2752–2762 (2005). 3. R. Lindstrot, R. Preusker, and J. Fischer, “Empirical correction of stray light within the MERIS oxygen A-band channel,” J. Atmos. Ocean. Technol. 27(7), 1185–1194 (2010). 4. M. di Bisceglie, R. Episcopo, C. Galdi, and S. L. Ullo, “Destriping MODIS data using overlapping field-of-view method,” IEEE Trans. Geosci. Rem. Sens. 47(2), 637–651 (2009). 5. M. Bouali and S. Ladjal, “Toward optimal destriping of modis data using a unidirectional variational model,” IEEE Trans. Geosci. Rem. Sens. 49(8), 2924–2935 (2011). 6. J. H. Ryu, H. J. Han, S. Cho, Y. J. Park, and Y. H. Ahn, “Overview of geostationary ocean color imager (GOCI) and GOCI data processing system (GDPS),” Ocean Sci. J. 47(3), 223–233 (2012). 7. W. Kim, J.-H. Ahn, and Y. J. Park, “Correction of stray-light-driven interslot radiometric discrepancy (ISRD) present in radiometric products of geostationary ocean color imager (GOCI),” IEEE Trans. Geosci. Rem. Sens. 53(10), 1–15 (2015). 8. K. E. Moore, M. G. Nicholson, and J. R. Warwick, “User-defined objects for stray-light and illumination system modelling,” Proc. SPIE 5249, 400–407 (2004). #261057 Received 15 Mar 2016; revised 19 Apr 2016; accepted 24 Apr 2016; published 2 May 2016 © 2016 OSA 16 May 2016 | Vol. 24, No. 10 | DOI:10.1364/OE.24.010232 | OPTICS EXPRESS 10232 9. G. L. Peterson, “Stray light calculation methods with optical ray trace software,” Proc. SPIE 3780, 132–137 (1999). 10. J. O. Park, W. K. Jang, S. H. Kim, H. S. Jang, and S. H. Lee, “Stray light analysis of high resolution camera for a low-earth-orbit satellite,” J. Opt. Soc. Korea 15(1), 52–55 (2011). 11. S. Seong, J. Yu, D. Ryu, J. Hong, J. Y. Yoon, S. W. Kim, J. H. Lee, and M. J. Shin, “Imaging and radiometric performance simulation for a new high-performance dual-band airborne reconnaissance camera,” Proc. SPIE 7307, 730705 (2009). 12. S. Jeong, Y. Jeong, D. Ryu, S. Kim, S. Cho, J. Hong, S. W. Kim, and H. S. Youn, “In-orbit imaging and radiometric performance prediction for flight model Geostationary Ocean Color Imager,” Proc. SPIE 7452, 74520F (2009). 13. E. Oh, S. W. Kim, Y. Jeong, S. Jeong, D. Ryu, S. Cho, J. H. Ryu, and Y. H. Ahn, “In-orbit image performance simulation for GOCI from integrated ray tracing computation,” Proc. SPIE 8175, 81751H (2011). 14. D. Ryu, S. Seong, J. M. Lee, J. Hong, S. Jeong, Y. Jeong, and S. W. Kim, “Integrated ray tracing simulation of spectral bio-signatures from full 3-D earth model,” Proc. SPIE 7441, 74410A (2009). 15. D. Ryu, S. W. Kim, D. W. Kim, J. M. Lee, H. Lee, W. H. Park, S. Seong, and S. J. Ham, “Integrated ray tracing simulation of annual variation of spectral bio-signatures from cloud free 3D optical Earth model,” Proc. SPIE 7819, 78190E (2010). 16. D. Ryu, S. W. Kim, and S. Seong, “Improved atmospheric 3D BSDF model in earthlike exoplanet using raytracing based method,” Proc. SPIE 8521, 85210F (2012). 17. D. Ryu, S. W. Kim, and S. Seong, “Ray tracing simulation of spectral bio-signatures with integrated earth BRDF model,” Proc. EPSC-DPS Joint Meeting, 631 (2011). 18. J. Yu, D. O. Ryu, S. H. Ahn, and S. W. Kim, “Real scale ray-tracing simulation of space earthshine measurement with improved BRDF model of lunar surface,” Proc. SPIE 8146, 81460X (2011). 19. I. Breault Research Organization, Sources in ASAP (Breault Research Organization, Inc., 2008). 20. “Standard solar constant and zero air mass solar spectral irradiation tables (ASTM E490 00a)”, retrieved http://www.astm.org/Standards/E490.htm. 21. D. C. Look; J. Look, “Diffuse reflection from a plane surface,” J. Opt. Soc. Am. 55(12), 1628–1632 (1965).


Introduction
In practical fields of remote sensing science, the suppression of stray light is required for image correction. Stray light from an unknown source may affect radiometric accuracy. Although stray light sources are considered before remote sensing instruments are launched, the unexpected presence of such light often impacts the acquired images. Therefore, when stray light phenomena have been discovered on remote sensing images, attempts have been made to apply corrections. For example, in the case of the FY-2 satellite, stray light is introduced from outside of the Earth's atmosphere and from direct solar irradiation at midnight. An analysis of the mechanism and modeling of direct reflection was performed, and a mathematical model based on images acquired daily was developed to reduce the impact of stray light on the FY-2 satellite [1]. The Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) is also affected by stray light from outside sources or the instrument itself, and this impact has been analyzed using calibrated in-orbit lunar images [2]. For another spaceborne spectrometer, the Medium Resolution Imaging Spectrometer (MERIS), an empirical correction method is used to reduce the stray-light radiometric errors caused by multiple reflections and scattering [3]. For the Moderate Resolution Imaging Spectrometer (MODIS), radiometric discrepancy phenomena such as striping noise, Sun glint, and stray light effects have been identified and studied. Striping noise is induced randomly by poor gain-and-offset calibration among detectors, scanning mirror operation error, and brightness effects in the observation target related to high reflectance, which occurs at levels close to sensor saturation [4,5]. The Geostationary Ocean Color Imager (GOCI), which was a target optical system in this study, has also been affected by a radiometric discrepancy issue due to stray light. The GOCI has a three-mirror anastigmat optical system [ Fig. 1(a)] with a 140-mm entrance pupil diameter. The filter wheel for a total of eight spectral bands is located in front of the detector, and a wedged filter system, which pairs of filter are combined with a different gap angle for each spectral band, is used to reduce ghost phenomena. The detector module features a complementary metal-oxide-semiconductor array of about 2 megapixels, with a pixel size of 14.81 μm (east-west) × 11.53 μm (north-south) to compensate for the Earth's oblique projection over the Korean Peninsula. The 2500 × 2500-km coverage area of the GOCI is shown in Fig. 1(b), with the sequence of slot acquisition for 16 slots centered at 130°E, 36°N, acquired sequentially using a step and stare observation method. Each slot has an overlap area of about 80 pixels and an acquisition time of about 2 minutes per image, which is due to the time required to operate the pointing mirror (POM) and filter wheel mechanism (FWM) when taking images at different bands [6]. Although the GOCI optics were designed to minimize the ghost effect, a radiometric issue referred to as the inter-slot radiometric discrepancy (ISRD) is present due to the effects of stray light, as shown in Fig. 2 [7]. Figure 2 shows the ISRD phenomena at the boundary between slots 7 and 10. The red arrows in images from bands 6 and 8 indicate the distinct lines separating the upper and lower slots, with a difference in radiance. Although the effect of radiance variations due to the change in the Sun's elevation angle is removed, the ISRD phenomenon is less than 3% in most overlap areas, but may be up to 15% for particular bands [7]. These conventional studies were performed with the aim of improving remotely sensed images by correcting the stray light effect based on mathematical and experimental methods, rather than identifying the reason for the effect. On the other hand, conventional stray light analysis conducted during the development phase is usually based on ray tracing, with consideration of the optical and mechanical geometry of the instrumental system. Optical properties such as reflectance, transmittance, and the bidirectional reflectance distribution function (BRDF) have been used to determine abnormal ray paths [8][9][10]. Furthermore, after launching and operating the satellite, it is quite complicate to simulate the similar operating environments with only controlling the illumination ray source. In short, ray trace-based stray light analysis has concentrated on the internal effects of optical systems and the ideal design of models; the estimation of in-orbit performance and the amount of stray light that could affect an image is thus difficult.
Previous studies of stray light issues using the aforementioned techniques have focused on optical systems per se, without consideration of unexpected external factors such as the operating mechanism acquisition concept, observational environment, and so on. To achieve not only the optimal design engineering, but also radiometric image processing, an in-orbit stray light reduction technique is needed, which could be applied to remote sensing research and used to understand the cause and effect of stray light on images. In this study, a new stray light analysis concept that considers the observational environment was developed. The method was based on an integrated ray tracing (IRT) technique, by which "ray tracing" was used for the simultaneous computation of imaging and radiative transfer equations. The IRT technique is used widely for end-to-end performance simulations for satellite instruments [11][12][13], Earth spectral analysis modeling [14][15][16][17], and Earthshine modeling between the Earth and Moon [18]. To embody the operational concept of the GOCI, we constructed an IRT model consisting of a Sun model as a light source, a target Earth model, and the GOCI optical system model. We then combined these models in a Monte Carlo-based ray tracing computation.
The remainder of this paper is organized as follows. Section 2 describes the IRT model for stray light analysis, including relevant concepts and a simulation procedure. Section 3 shows the results of simulation using the IRT method for the GOCI, based on two case studies. ISRD phenomena candidates for the GOCI images and the limits of the IRT method are discussed in section 4. Finally, section 5 provides a summary of the new method and the conclusions of the study. The new ray tracing technique was intended to combine the light source (phase 1), measurement target (phase 2), and instrument structure (phase 3) in a single Monte-Carlo ray traction computation procedure, as illustrated in Fig. 3(a). We saved the ray tracing information for the flux and direction in each step, and then generated a new ray set in the next step [ Fig. 3(b)]. The source model characteristics were changed to account for light emitted from the Sun (source) and traveling to Earth (target), and reflection and scattering of the rays at the surface of the Earth (target) to the instrument in orbit. The following three equations were used to obtain ray information for the simulation in each step. (

A new concept of stray light analysis based on IRT
In phase 1, rays from the Sun were generated and emitted in all directions in a statistical manner, and then traced to the target area on Earth. The rays were then scattered and reflected, and the intensity ( ( ) In phase 1 of this study, the initial source model for Sun scattering was determined with a random polar grid source, generating the calculated power for the particular simulation circumstance [19]. We assumed that the Sun is a Lambertian spherical source, with a total energy of 591.4 W/m 2 at a wavelength range of 400-900 nm (defined as the solar constant) [20]. In the Sun model construction, the solar radiance L s and solar radiant flux Φ S at the Sun's surface were computed from the solar irradiance data E S(earth) at the top of the atmosphere (TOA) using Eqs. (4) and (5).
where D SE ( = 1.4960 × 10 11 m) is the Sun-Earth distance and R S ( = 6.9550 × 10 8 m) is the Sun's radius. The flux reaching the Sun's observation target area was calculated using Eq. (6): where E S(TA) indicates the solar irradiance data at the target area (A S(TA) , 6.2500 × 10 12 m 2 ), which was the area observed by the GOCI. The reflected intensity I E(TA) was calculated using Eq. (7). The equation was derived from the general expression for Lambert diffuse reflection [21].
where , TOA E λ ρ is the reflectivity of each surface property, R E ( = 6.3710 × 10 6 m) is the Earth's mean radius, φ is the azimuth angle, ε is the angle between the Sun and the GOCI azimuth angle, θ is the latitude angle, and α is the angle between the Sun and the GOCI latitude angle [ Fig. 4.]. To calculate ( ) E TA I in this study, we adjusted the three simulation cases using the orbit coordinate parameters, as summarized in Table 1. The azimuth and latitude angles of the Sun (θ SUN , φ SUN ) were derived from the ephemeris values according to the time and day, which was the same as in the real GOCI image acquisition environment, and we used standard orbit data from the GOCI (θ GOCI , φ GOCI ) because it is a geostationary satellite located in the same position. In addition, we assumed the TOA reflectance ( , TOA E λ ρ ) over the sea surface to be 0.7, and land and cloud TOA reflectance to be 0.3 and 0.99, respectively, in the IRT simulation. ), were assumed that the rays traveled in all spatial directions, with optical properties determined by the instrument. During this process, we had to control the ray information (intensity amount and ray direction vectors) heading in the instrument's direction. Then, Φ E(TA) (W) and Φ G(AP) (W) were calculated using Eq. (8) and Eq. 9: where D EG is the Earth-GOCI distance ( = 36,000 km), which was fixed in this study. The relative contribution of each path's flux to all returned rays. c The rear surface faces the detector along the line of sight, and the front surface faces the opposite way.

Phase 3: in-depth stray light analysis at the instrument level ( DET Φ )
Finally, in phase 3, to estimate the flux on the detector plane ( DET Φ ), we completed two processes. The first process was the identification of stray light and ghost candidates in the instrument and adjustment of each optical property to the candidates. The second process was to create a simulation with the same operational activity of the instrument. In the case of the GOCI, we followed the same operational procedure and controlled the pointing mirror angle to acquire the flux on the detector plane, with the simulation image. Table 2 shows the ray paths of the illuminated and critical objects which were suspected to be the candidates of stray light sources. Although the GOCI optical system lacks a lens element, each pair of spectral filter in the filter wheel can generate split rays by internal reflection and make a ghost image at the detector plane. The locations where the ray paths reached a thermal heater and a support beam are illustrated in Table 2. The path ratio reflected from the wide and planar thermal heater, which faced parallel to the focal plane, covered more than 93% of the critical path and 94% of the illuminated path. In the case of a support beam, the scattered rays from M3 or M4 reached about 7% of the critical and illuminated paths.
Multi-reflected rays between the front and back surfaces of filters, another stray light source, were simulated at the instrument level. In the simulation, we focused on the signal impacts from light sources outside of the field of view (FOV). Generally, one GOCI slot image covers ± 0.4° in the north-south direction. In the simulation, we set the boundaries of FOV coverage at 0.4° to 1.2°, and executed each ghost analysis according to an incident angle of 0.03°. Figure 5 shows the ghost signals (excluding nominal signals) on the detector plane at each incident angle. Figure 5(a) shows that no ghost signal from sources outside of the FOV was present until the incident angle was 0.7°. Most ghost effects were revealed in the incident angle range of 0.78-0.90°, as shown in Figs. 5(b)-5(f). The light sources outside of the upper slot area of the FOV affected the detector plane in the southward direction. These results indicate that a bright light source located in the center of the upper slot (about 0.8-0.9°) would influence the lower area, below the slot image, in the southerly direction. These instrument-level simulations were based on an in-depth IRT analysis simulation, which is described in the following section. Following the real operational sequences of the GOCI, a simulation was conducted in which the satellite was targeting a 5.1° area. Each subsequent simulated slot image was controlled using the POM angle on the GOCI instrument. To simulate with the same operational environment, we calculated the simulation parameters (nominal vector angle on each axis) of POM angles with regard to the slot center positions on the GOCI image, as shown in Table 3 [13].
Considering the complex pointing mirror mechanism of the step-and-stare observation technique used for the GOCI [6], the IRT technique has the advantage of calculating the relationship between the center of the target area and the rotation angle of the pointing mirror at a specific time, with the specific Sun altitude relative to the target area and the instrument orientation. The IRT method described here can be used to build a well-composed simulation condition that is similar to the real in-orbit environment. Furthermore, to determine how unwanted light affects the acquired image, the IRT technique is useful for examining the amount and origin of stray light. The verification results are shown in the following section.

Case 1: stray light pattern in a cloud-free environment
For the three simulation phases, the first case considered the effect of stray light under cloudfree conditions. We compared simulated images according to the acquired time variation in TOA-level radiance. Figure 6 shows the stray light pattern that was subtracted from a mosaicked image by selecting abnormal ray paths in spectral cases. We selected and compared simulation images among bands 2, 6, and 8. Generally, these stray light patterns show the degradation pattern in the north-south direction in each slot image. The stray light pattern changed according to the time variation, as shown in Fig. 6. The UTC-00 simulation results showed that the effect was greater in the southeastern area than in other locations. On the other hand, the UTC-06 results indicated that a bright stray light pattern occurred in the southwestern area. The lowest slot areas (slots 13-16) had a wide zone affected by stray light, which reached about 30%, as shown in Fig. 6(e). Table 4 shows mean, minimum, and maximum periodical and spectral stray light ratios. The minimum values were used to match the nominal background effects of stray light, and the maximum values indicate the strongest stray light effect in slot areas, which were usually located in the southern areas. The mean values, including the standard deviations, indicate the representative stray light effects in the total slot area (i.e., the area of slot 10 in Table 4). The largest amount of stray light at UTC-03 was 1.3340% at band 2 (443 nm) in the area of slot 10 via the nominal radiance in the cloud-free simulation at all spectral bands. Spectral comparison showed that the simulation images from near-UV bands (bands 1-3) had relatively large background stray light effects. Standard deviations for bands 5 and 6 were small, indicating that the region affected by stray light was spread over most slot areas.   Figure 7 shows the variation in the effect of stray light from a bright source outside of the FOV. When the cloud source was located in slot 8, the most powerful effect of stray light was in the same location as the cloud, but the semicircular stray light pattern appeared at the slot 9 position. These results indicate that the southern slot can be affected from a northern bright source as a ghost phenomenon. In addition, a cloud-induced ghost pattern followed a cloud moving from west to east. The spectral stray light ratios induced from a nearby bright source are summarized in Table 5. The minimum values are almost the same as those shown in Table 4, but the mean and maximum values of all spectral images were higher due to the brightness of the source. In the cloud-free simulation, the largest effect was observed for the UTC-03 results; in the cloudy simulation, the greatest effect was found for the UTC-06 results. The maximum stray light effect was 8.5841% in the band 2 image at UTC-06. This amount was 6.65 times greater than that of the cloud-free simulation result. Figure 8 shows stray light profiles for different bands in cloud-free and cloudy simulations under the same conditions (slot 10 area, UTC-03). Each colored dot is an average of the 10 surrounding pixels, and the error bar indicates its standard deviation. In the cloud-free simulation, an abnormal stray light pattern was apparent from about 3600 lines, and those values gradually increased. On the other hand, we noticed that the cloud-induced stray light effect region started at about 3800 lines, in which the pixels were about 1/3 of the distance from the southern area of slot 10.    [16]. (c, f) Stray light-corrected images obtained by IRT simulation.

Discussion
As already mentioned in the previous section, simulation results show that stray light appear at the lower end of slot boundary in a clear sky condition [Fig. 6], and a bright light source located on upper slot also makes stray light as a ghost effect [ Fig. 7]. Although we connot sure that it is the only reason for the present ISRD problem in the GOCI data, We attempted to apply the simulated results, obtained with the proposed method, to the results of other studies for verification [7]. Figure 9 shows two spectral band images from 680 and 865 nm [Figs. 9(a) and 9(d)], which are compared with the stray light effect determined by the method proposed by Kim et al. (2015). In the real acquired images [Figs. 9(a) and 9(d)], we could figure out the ISRD pattern with a complicated gradation shape in about 400 pixel columns at the lower end image [ Fig. 9(b) and 9(e)]. This ISRD pattern is expected that diffused stray light bias made the largest contribution to the differences among slot boundaries, which is similar to that seen on simulation results. By examining correspondence between the simulated stray light pattern and that appearing on real acquired images, we removed the stray light effect at each pixel value with based on simulation results, and then obtained first-order corrected images, as shown in Figs. 9(c) and 9(f). Although not all effects of stray light were removed, they were drastically decreased in the bottom slot areas on both images. For example, the red patterns in Figs. 9(b) and 9(e) changed to green, indicating the absence of a stray light effect. The stray light ratio in Fig. 9(b) has a maximum value of 2.54% and mean value of 2.38% ± 0.16% in the dotted box area. After correction using the IRT method, the stray light ratio was reduced to 1.02% ± 0.17% and the maximum value was 1.21% [ Fig. 9(c)]. The results for the band 8 case also showed that a decrease in the average value from 1.09% ± 0.07% to 0.35% ± 0.07%, and a decrease in the maximum value from 1.18% to 0.47%, as shown in Figs. 9(e) and 9(f). This verification results indicate that IRT based simulation have generated statistically similar pattern comparing with ISRD phenomena in a real acquisition image. Compared with conventional stray light analysis methods, the advantages of the proposed method are as follows: 1) it can be applied to most remote sensing environments after payloads have been successfully launched in orbit; 2) as a ray tracing-based method, it can easily determine stray light candidates; and 3) the inner (instrument-level) and outer (observational environment-level) conditions can be modified for in-depth applications. Its limitation is that the current analytical environment is focused on TOA-level radiance, with no consideration of atmospheric conditions. Thus, quantitative comparison with real acquired remote sensing images is difficult, and direct correction is currently challenging.

Conclusion
In this study, we developed a new stray light analysis method for step-and-stare type instruments in geostationary orbit that incorporates a Monte Carlo based ray tracing technique, which can be used to determine historical and radiometric reasons for stray light.
The results proved that a changeable bright source was a candidate source of stray light and that it disturbed the closure of the slot image from the next slot target. The IRT end-to-end simulation technique revealed evidence of ISRD phenomena during the use of the GOCI. Elements of the observational environment, such as the orbit, Sun's zenith angle, and instrument characteristics, were fully considered to determine the stray light path. The IRT technique developed in this study not only overcomes the limitations of these conventional stray light analysis methods by modifying the details of the Sun and Earth models and applying instrument parameters directly, but it also actively controls the orientation of the simulation environment and traces the ray paths of stray light candidates, which helps to determine the reason for stray light effects.
In conclusion, stray light analysis is usually used in the development phase to determine the optimal design of a sensor. However, the in-orbit stray light analysis method developed in this study confirmed that it is also useful for the identification of discrepancies or abnormal phenomena in remotely sensed images, and for supporting the image correction algorithm. The IRT method presented in this paper can be used to build well-composed simulation conditions similar to the real in-orbit environment. Furthermore, to determine how unwanted light affects an acquired image, the IRT technique can be used to determine the amount of stray light and its origin. In the near future, the specific stray light correction algorithm will be further investigated.