1. Introduction

The optical performance evaluation of space instruments is an important factor in successful scientific missions involving remote sensing of planetary objects such as the Earth and other solar planets, satellites, and even exoplanets. These evaluations have been accomplished in several ways, including physical experiments and end-to-end simulations during the instrument design phase. However, the lack of performance evaluation criteria has caused issues during the on-orbit operations of space missions. Examples include the spherical aberration of the Hubble Space Telescope (HST) [1] and imaging defocus of the High-Resolution Instrument (HRI) of Deep Impact/EPOXI [2]. These issues, regarding imaging optics in particular, have been identified by estimating the point spread function (PSF) using starlight observation or imaging simulations of a point source; i.e., Tiny Tim with the HST [3] and the research by Barry et al. for the EPOXI mission [4]. However, even though these methods were used to observe the Earth and other planets, they did not directly simulate the Earth or other planets to establish a measurement target.

Most typical planetary remote sensing missions, including those from the Earth to exoplanets, have models for estimating instrument performance during on-orbit operation. These models usually do not allow for other optical instrument models, except in the case of a corresponding instrument involved in a different mission. In addition, several types of physical instrument parameters are limited by the scientific parameters of the mission. For example, the NPOESS simulator [5] computes a radiance distribution but does not measure imaging or spectral performance. PICASSO [6], S2eteS for Sentinel-2 [7], and the model proposed by Bu et al. [8] calculate remote sensing images for instruments, but they have limitations associated with their ability to calculate spectral results. In contrast, SpecSim for JWST MIRI IFU [9], AxeSim for HST WFC3 [10], and ExoSim [11] calculate the transit spectroscopy of exoplanets but do not provide imaging results. Hyperspectral imager simulators (e.g., the simulator for CRISM [12], E2ES for FLEX FLORIS [13,14], EeteS for EnMAPs [15], and HITESEM [16]) can calculate both imaging and spectroscopy parameters, but some of these instruments are not compatible with other instruments.

In simulating the Earth or other planets as measurement targets, some instrument models utilize radiative transfer atmospheric models and surface reflectance models, i.e., the simulators for NPOESS [5], CRISM [12], and DIRSIG [13]. These models include the local area of the planet, which is limited by the field of view (FOV) of the instrument in use, rather than the global planetary surface. The DIRSIG [16] model has the most powerful scene generator, which imports three-dimensional (3-D) freeform objects as measurement targets. However, it is limited in its ability to evaluate the imaging performance of an optical instrument alone. Others (i.e., EeteS [15], S2eteS [7], PICASSO [6] and Bu’s model [8]), which do not have a model measurement target, require highly resolved images from other observations. These models utilize a PSF convolution using a high-resolution image and an analytical atmospheric correction. If the input data are of lower resolution than the instrument’s output data or if there are no input data, as in the case of undiscovered exoplanets, it is difficult to evaluate the performance of the optical instrument.

Global earth system models have been developed to perform spectroscopy on Earth-like exoplanets. Early earth system models use the optical characteristics of the Earth, such as its absorption, transmission, reflection, refraction, and scattering effect. For example, Ford et al.’s model lacks an atmosphere or cloud model that interacts vertically with the surface [17]. Montanes-Rodriguez et al.’s line-by-line radiative transfer model only averages the one-dimensional atmospheric profile horizontally; it does not calculate spatially resolved images or spectra [18]. The three-dimensional line-by-line or band model is designed for some radiative transfer codes of terrestrial atmospheres [19–21]. This algorithm leads to a 3-D earth system model composed of time and disk-resolved input data from remote sensing observations. Some recently developed earth system models calculate 3-D line-by-line spectral results for the Earth with its optical characteristics [22–26]. However, these earth system models do not consider instrument imaging performance for the Earth.

An integrated ray tracing (IRT) model was proposed in 2002 [27]. This model uses a global earth climate model with a uniform Lambertian scatter surface and an Amon-RA imaging instrument for measurement of the albedo of the Earth. The first earth system model to have six types of Lambertian scatter surfaces and a clear-sky atmosphere, including specular transmission, reflection, and absorption except for scattering, was developed in 2009 [28]. Since 2010, non-Lambertian scatter models have been imported into clear-sky earth system models, i.e., the single scattering of the atmosphere, the sun glint scattering of the ocean's surface, and the parametric bidirectional reflectance distribution functions (BRDF) of six land surface types [29–31]. Previously developed earth system models that do not include a cloud model are limited to use in comparing simulation results with corresponding observations for clear-sky conditions only. Therefore, the performance evaluation of space optical instruments is limited in the case of cloudy-sky atmospheric conditions on the Earth.

The limitations of the existing models discussed above highlight the need for a new, combined global earth system and instrument model for use in evaluating the optical performance of remote-sensing space instruments such as imagers, spectrometers, and hyperspectral imagers. Such a model should be able to calculate imaging and spectroscopic results simultaneously. Moreover, a global earth system model used as a scene generator for measurement targeting should provide fully resolved spatial, spectral, and temporal coverage of the Earth, with sub-models for the molecular atmosphere, clouds, aerosols, land surfaces, snow cover, sea ice, and the ocean’s surface. In particular, a cloud model that can be used for performance evaluation in cloudy sky conditions should be included in the new earth system model.

The conceptual basis of the IRT method and the IRT ray tracing algorithm are introduced in Section 2 of this paper. In Section 3, the new global earth system model is described and is shown to improve on not only the three types of cloud sub-models but also differentiation of the 16 types of land surfaces, using a semi-empirical BRDF based on remote sensing measurements a spatially and temporally resolved input database. Section 3 also describes a model for a diffraction-limited instrument with variable specifications for different levels of optical performance. In Section 4, imaging simulation results obtained with the new earth system model are compared with EPOXI HRI-VIS measurement data. The adjustable-performance instrument model was used to evaluate the optical performance of HRI-VIS using a direct comparison between the imaging simulation results and the observations. Finally, a line-by-line spectral simulation was used to derive disk-averaged spectra that were compared with other spectra from the traditional radiative transfer model used in previous studies. Discussion of the study results and conclusion that may be drawn from the results are presented in Sections 5 and 6, respectively.

2. Conceptual basis of integrated ray tracing (IRT) method

2.1 Ray tracing method

The basic elements of the IRT method are the rays, media, and objects that are represented in Fig. 1. A ray is defined via a ray tracing algorithm, according to Eq. (1) [32]:

 

Fig. 1 Illustration of specular (A) and scattering (B) ray tracing.

Download Full Size | PDF

As Fig. 1(a) shows, monochromatic rays are emitted in empty space (object number o = 0) with an initial medium (m = 0). Then, traced rays reach and stop at the first generated object (o = 1) in the case of an absorbing sphere. If the first object is on the border of another medium (m = 1), the rays are traced in the specular direction:

As Fig. 1(b) shows, if the object surface is a scattering layer and the extra j rays are emitted in the direction of reflection and transmission,

For efficiency in the ray tracing calculations, an importance edge surface is defined. The position (${\overrightarrow{P}}_{S,j}$) and direction (${\overrightarrow{D}}_{S,j}$) of the new rays are randomly defined at a selected solid angle from the importance edge surface. In the transmitted scattering case, rays are refracted out from the importance edge based on Snell’s law.

Scattering radiant power, F_S,j, and integrated scattering power for all directions (S) are calculated using the following equations [32]:

This concept of ray tracing is implemented in the Advanced System Analysis Program^TM (ASAP) ray tracing optical simulation software available from the Breault Research Organization [32]. This industry-standard software is widely used in optical system analysis. The functionality of this software was expanded in this study for use as an end-to-end simulation tool using the proposed earth system model.

2.2 Flow of the integrated ray tracing (IRT) method

The IRT method was developed for use in an “end-to-end” simulation, where the first “end” term is a light source for illuminating a measurement target and the other “end” term is the output data of the measurement target from an instrument detector with optical performance [33].

When running end-to-end simulations using the IRT method, the input data from an input database is assigned to sub-models in the earth system and instrument models, as shown in the area enclosed by the gray dashed line in Fig. 2. The input database consists of the spatial, directional, spectral, and temporal information required for the sub-models and instrument model used in the IRT method. Table 1 summarizes the input database’s information, including the data resolution, the numbers of media and objects, and the optical properties of the sub-models.

 

Fig. 2 Flowchart of the IRT method showing molecular atmosphere sub-model enveloping other sub-models within the earth system model.

Download Full Size | PDF

Table 1. Resolution of Input Data and Optical Properties of Objects in Sub-Models

View Table | View all tables in this article

The ray tracing process is shown as a red arrows flowchart in Fig. 2. The sub-model of the sun used as a light source emits and traces the rays toward the molecular atmosphere sub-model that encircles the other sub-models of the earth system model. The rays from the sun model are absorbed, refracted, transmitted, and reflected by molecular layers in the specular direction in the atmosphere system model. Next, the rays are traced to clouds, aerosols, land covers, ocean, snow cover, and sea ice sub-models through the molecular atmospheric layers. Then, all of the sub-models reflect and scatter the rays toward the entrance pupil of the instrument's aperture, thus leaving the molecular atmosphere sub-model.

In this study, the instrument model used was a hypothetical Cassegrain telescope imager with diffraction-limited performance. The rays at the entrance pupil were traced to the optical parts in the instrument models. The non-sequential tracing was terminated as the traced rays were absorbed in the focal plane or by opto-mechanical structures. The image and spectrum of the earth system model can be calculated simultaneously using the position, direction, and radiant power of the rays at the focal plane after ray tracing is performed and the data are saved in an output database. The traced rays at the top of the atmosphere with regard to the sun and the entrance pupil of the instrument are also saved in the output database.

3. Details of sub-models

3.1 Sun

The spectrum of the sun sub-model was produced using Gueymard’s solar spectral irradiance model [34] and 6-h-scale total solar irradiance (TSI) observation data from the Solar Radiation and Climate Experiment (SORCE) mission [35]. The TSI of Gueymard’s spectrum was based on a solar constant of 1366.11 W/m² [34], whereas that of the SORCE mission observations ranged from 1360.5256 to 1360.4699 W/m² for the simulations [35]. Therefore, a radiometric correction was performed to match the TSI of the spectrum model to that of and observations. A total of 30,000 lines of spectral irradiance were interpolated to 295 lines with a wavelength range from 300 nm to 2500 nm. Additionally, four temporal data points of 6-h-scale TSI were interpolated to 25 segments with a 1-h scale.

The spectral and temporally resolved solar irradiance (W/m²/nm) at a distance of 1 AU was converted to a hemispheric spectral solar radiant power (W/nm). This value was set for emitted rays in the sun sub-model hemisphere, with a radius 6.957 × 10⁸ m [36]. One hundred thousand rays were emitted per unit wavelength and time for the Earth disk-averaged spectra simulation. The imaging simulations used 9,000,000 input rays to reduce ray tracing noise. The surface of the sun sub-model was assumed to have a Lambertian scattering nature, according to which rays scatter toward the importance edge of the earth sub-models, 0.995572 AU away [37] from the sun sub-model under the initial conditions.

3.2 Molecular atmosphere

The molecular atmospheric sub-model consists of two parts, as shown in Fig. 3: (1) atmospheric profile production and (2) modeling of atmospheric optical properties. A Tape5 input file for the Line-By-Line Radiative Transfer Model (LBLRTM) is generated by the Tape5 generator, which is the interface for running the LBLRTM, by importing the initial pressure (P), temperature (T), and gas mixing ratio (mr) profiles, height (H) selection data, and wavelength boundaries [38]. The heights of the layers were selected to be 0.4, 0.8, 1.0, 1.5, 2.0, 3.0, 3.5, 5.0, 6.0, 10, 12, 25, 35, 50, and 100 km above the WGS84 [39] ellipsoid. The pressure profile, temperature profile, and gas mixing ratio data were obtained from the LBLRTM default atmosphere.

 

Fig. 3 Flowchart of the molecular atmospheric sub-model.

Download Full Size | PDF

Next, the LBLRTM loads the Tape5 input condition file and Tape3 line database file and calculates the spectral absorption and scattering coefficient profiles. The optical depth profiles of the atmosphere are calculated using the coefficient profiles multiplied by the optical path length. The optical depth is used as the molecular extinction, using Beer’s law to attenuate the radiant power of rays passing through the atmospheric layer. The LBLRTM calculates both coefficient profiles and auxiliary output profiles; i.e., pressure, temperature, and air density profiles. The spectral refractive indexes of the atmospheric layers are calculated from the pressure, temperature, and air density profiles using the LOWTRAN 7 equation [40]. The rays at the atmospheric layers are refracted based on their refractive indexes, according to Snell’s law. While being refracted and scattered, the directions of escaping rays from the outermost layer are adjusted to the instrument aperture using an additional importance edge technique [32]. The Fresnel equation is used to calculate the specular reflection of the atmosphere using refractive indexes. The transmission of specular rays was calculated by subtracting the molecular extinction and specular reflection from Fresnel reflection from the radiant power of the incident rays. The scattering model for the molecular atmosphere is based on Rayleigh scattering and is derived from the “successive order of scattering” method, represented by the following equations [41].

Finally, the molecular atmospheric optical model combines the calculated extinction, specular reflection, specular transmission, and Rayleigh scattering of the molecular atmosphere. This model is based on the imported molecular atmospheric structure, which is constructed with 15 vertical layers and five latitudinal distributions from the LBLRTM's default five latitudinal models—tropical, arctic summer, arctic winter, mid-latitude summer, and mid-latitude winter [20,38]—on a WGS84 ellipsoid.

3.3 Clouds and aerosols

The three layers of the cloud sub-model, including the low, middle, and high clouds, were assigned to the third (1.0 km), eighth (5.0 km), and tenth (10 km) layers of the atmosphere, respectively [42]. The cloud optical model calculation process is represented in Fig. 4. The liquid water content (LWC) [43] and effective radii [43] of the water particles were used as inputs to the parametric model for a liquid water cloud proposed by Hu and Stamnes [44]. This model was used to calculate the scattering coefficient, single scattering albedo, and asymmetric factor of the liquid water cloud for the low and middle levels. The database for the ice water cloud model proposed by Fu et al. [45] was used to calculate the scattering parameters of the ice water cloud for the high-level cloud, using the ice water content (IWC) [45] and generalized effective size [45].

 

Fig. 4 Flowchart of cloud sub-model.

Download Full Size | PDF

Next, the single scattering cloud phase function was evaluated by substituting an asymmetric factor into the Henry–Greenstein (HG) phase function. Then, the effective phase function of the cloud was derived using the equation proposed by Pfeiffer and Chapman for multiple scattering estimation [46,47], which estimates the multiple scattering of clouds. Next, the single-layer optical model, including multiple scattering by clouds, was constructed with the effective phase function, scattering coefficient, and single scattering albedo by the “successive order of scattering” method [41]. This optical model was used in the ten-level discrete fraction map of the clouds, which was constructed from the low, middle, and high cloud fraction data in the MERRA-2 retrieval database for a 1-h timescale [48]. The ten-level fractions were separated by 0.995, 0.975, 0.925, 0.850, 0.750, 0.600, 0.400, 0.200, 0.050, and 1E-8.

Five types of aerosol sub-models, including continental average, maritime clean, desert, mineral transported, and stratosphere models, were assigned to the fifth (2.0 km), seventh (3.5 km), and thirteenth (35 km) layers of the atmosphere model [49]. Three types of aerosols (continental, maritime, and desert) were modeled in a single layer at 2.0 km and were mapped with MODIS land cover data. The aerosol optical model calculation process is similar to that for the cloud model, but the aerosol model uses an OPAC database [50] that provides the extinction coefficient, single scattering albedo, and asymmetric factor for each aerosol model.

3.4 Land, ocean, snow cover and sea ice

The calculation process of land, ocean, snow cover and sea ice sub-models is represented in Fig. 5. Using the “BRDF-albedo model parameters 16-day L3 0.05 deg CMG (MCD43C1)” in the MODIS product [51], the K₀, K₁, and K₂ data of seven bands (blue, 470 nm; green 555 nm; red, 645 nm; near-infrared (NIR), 860 nm; and three IR bands, 1.2 µm, 1.6 µm, and 2.1 µm [52]) were obtained from March 13 to 28, 2008 and were classified into the 16 types of land surfaces listed in “Land cover types yearly L3 global 0.05 deg CMG (MCD12C1)” [53]. Then, the seven-band BRDF data were spectrally interpolated to obtain line-by-line BRDF data (K₀(λ), K₁(λ) and K₂(λ)) using the spectra of conifers, deciduous trees, grasses, dry grasses, soil, wet soil, dune sand, and snow in the SMARTS radiative transfer code [54]. For some land surface types, two spectra were combined with different weights, as shown in Table 2.

 

Fig. 5 Flowchart of land, ocean, snow cover, and sea ice sub-models.

Download Full Size | PDF

Table 2. Land Surface Composition

View Table | View all tables in this article

For example, a land covered with woody savanna uses a combined spectrum of 60% grass and 40% soil in the spectral database of the SMARTS radiative transfer code. The combined spectrum is then scaled to fit approximately to seven-band BRDF measurements of the 8th land cover type. The fitted spectrum is converted to spectral K0, K1, and K2. Finally, these values are linearly interpolated to 295 lines with a wavelength range from 300 nm to 2500 nm.

The spectral BRDF of the land sub-model including ice and snow was defined by a semi-empirical parametric kernel model used for MODIS [55–57], expressed by Eq. (9). The parameters of this equation have been described in previous studies [29–31].

Land model structures were constructed by a Cartesian coordinate conversion from polar coordinate (longitude, latitude, and altitude) data on a WGS84 ellipsoid [39]. Longitude and latitude data from the MCD12C1 land cover map were used for 16 land surface types with a 1-yr timescale. The mean altitude of each surface type shown in Table 2 was calculated using digital topographic data from the “ETOPO1 Global Relief Model” published by the National Geophysical Data Center (NGDC) [58]. The monthly snow cover map from March 2008 from “Snow cover monthly L3 global 0.05 deg CMG (MOD10CM)” [59], and the daily sea ice map from March 18 to 19, 2008 from “AMSR-E/Aqua daily L3 12.5 km sea ice concentration” [60] were used to construct the snow cover and sea ice surface for the ocean sub-model and land sub-model. These two sub-models had the same snow BRDF as the land sub-model.

The ocean sub-model is modeled using Eq. (10) [29]:

The $BRD{F}_L(\lambda )$ term is defined by the ocean spectral albedo after subtracting the total integrated scattering of the sun glint scatter model ($BRD{F}_{SG}({\theta }_s,{\theta }_v,{\varphi }_s,{\varphi }_v,\lambda )$), according to Nakajima et al. [61], as described in previous studies [29–31].

3.5 Instrument model for imaging performance

In this study, we developed a hypothetical imaging instrument model for imaging purposes. The optimized imaging optics satisfied the diffraction-limited specifications of HRI-VIS data from the EPOXI mission [2]. The two-mirror Cassegrain telescope had a 300-mm-diameter aperture, a secondary mirror with spiders obscuring rays, and a 21.504 mm × 21.504 mm focal plane [2]. The telescope was designed to have a 10,500-mm effective focal length and a 0.118° FOV [2].

The root mean square (RMS) spot radii of the nine fields at the focal plane were approximately 1.7 to 5.9 μm. The spot sizes, modulation transfer function (MTF), and RMS wavefront errors (WFEs) were better than the diffraction-limited performance of geometric optics. The seven-band transmission reported by Hampton et al. [2] was applied to the filter window. The two mirrors were assumed to have 100% reflectance, with no scattering, no absorption, and no transmission from their surfaces.

Figure 6 shows a cross section of the instrument model used. The rays at the telescope’s aperture from the earth system model were traced to the focal plane via the optical parts of the instrument model. The telescope tube and instrument cover enclose the optical elements to block all rays except those in the field-of-view (FOV) direction. There are no internal baffles or vanes, but the traced rays that have abnormal paths due to scattering are excluded by selection using the tracing history. At the focal plane, quantum efficiency was adjusted for rays as a form of transmittance. Finally, an imaging result from the earth system model was obtained from the ray energy distribution. Neglecting filter transmission and quantum efficiency were also considered in comparing the simulated images with the calibrated images from HRI-VIS measurements obtained during the EPOXI mission. This comparison is shown in Fig. 7.

 

Fig. 6 Instrument design and traced rays.

Download Full Size | PDF

 

Fig. 7 RGB true-color image from simulation and observation. Two movie files from the simulation and one movie file from the observation over 24 h are included. Visualization 1 shows the simulation images of the Earth surface without atmosphere. Visualization 2 is the simulation for the Earth with atmospheric components. Visualization 3 shows the Earth images from the EPOXI HRI-VIS observation.

Download Full Size | PDF

In practice, the aberration in the performance of the HRI-VIS was detected using defocused imaging and was expected due to the uptake of humidity by the telescope before its launch [62]. Therefore, while controlling the distance between the primary and secondary mirrors, the defocused images were represented by a ray tracing simulation with the optical design that is described in Sections 4.3 and 5.

4. Case simulation and results

4.1 Case definition

In this study, the Earth-observing simulations were compared to EPOXI mission observations, including other research calculations, to verify the IRT method. The HRI-VIS of the EPOXI mission observed the Earth rotating for approximately 24 h between UTC March 18, 2008 at 18:19 and UTC March 19, 2008 at 18:19 [63]. A summary of this observation was presented by Livengood et al. [63], who termed it “EarthObs1.” The calibrated image data set of the Earth from the HRI-VIS data collected during the EPOXI mission was provided by the Small Bodies Node, which is a part of NASA’s planetary data system [37].

Three bands (350, 750, and 950 nm) from a total of seven wavelength filter bands in the HRI-VIS were acquired once per hour, and the other wavelengths (450, 550, 650, and 850 nm) were acquired every 15 min [37]. These acquired images were selected each hour, starting from the time the first 350-nm image was taken, for comparison with the simulation’s image.

The colors of the HRI-VIS filter were defined as the bands of violet at 350 nm, blue at 450 nm, green at 550 nm, orange at 650 nm, red at 750 nm, near-infrared at 850 nm, and infrared at 950 nm [2,37]. Typically, true-color RGB images have used the near-650-nm wavelength band as the red color. MODIS, for example, uses band 1 (620–670 nm) [64]. In this study, the orange band (650 nm) was selected for this purpose, instead of the red band (750 nm) in the HRI-VIS filter, to compose a true-color image. Table 3 summarizes the case details.

Table 3. Summary of Case Definition

View Table | View all tables in this article

4.2 Imaging result

Figures 7(a) and 7(c) show the true-color images, with a composite of orange, green, and blue bands, from the simulation obtained using the hypothetical HRI-VIS imager at two different times (18:19 on March 18, 2008 and 6:19 on the following day) [37]. The observation images from the HRI-VIS (Figs. 7(b) and 7(d)) [37] were used as references for verifying the imaging simulation. Similar to the imaging simulation, the single-band images are composed of the orange, green, and blue wavelength bands. Before the composition was performed, these single-band images were translated toward the center of the FOV to correct a pointing error of the HRI-VIS instrument.

The imaging simulation result was quite similar to the observation result. The simulation image was 178 pixels in equatorial diameter and 231 pixels in polar diameter, whereas the corresponding dimensions of the observation image were 194 and 243 pixels, respectively. These simulation measurements indicate that the surface scale was 55.0 km per pixel, whereas the scale of the HRI-VIS instrument is 54.9 km per pixel, according to Hampton et al [2]. This scale is almost the same value as that estimated from the simulation. The simulation and observation images were also similar with regard to the solar phase angle of approximately 57.7° and the 76% illumination fraction [63]. The distributions of the continent and ocean area in both images showed the same viewing geometry and rotation axis of the Earth. The colors of the ocean and land area regions that were not covered with clouds in the simulation images are comparable to those in observations. More detailed color results are described in the comparison of the spectral results presented in Section 4.4.

The sun glint scattering is shown via red circles on both the imaging simulation and observation images in Fig. 7. This scattering feature was disturbed by a thin cloud layer in the case of the simulation for March 19, 2008 at 6:18. Additionally, the distribution of bright clouds was similar, whereas the specific brightness distribution of the clouds, especially thin clouds, was different between the simulation and the observation. It is clear that the observation images are defocused due to instrument aberration in comparison to the simulation images. This phenomenon is discussed in detail in the next subsection.

4.3 Imaging performance evaluation

The defocused or blurred observation images were not compared directly to the diffraction-limited imaging simulation from the hypothetical HRI-VIS imager. However, as mentioned in Section 3.5, the instrument model used in the IRT method is able to modify the optical design, varying the distance between the primary and secondary mirrors to induce a defocus aberration. Furthermore, the optical design of the IRT instrument model is similar to the HRI-VIS instrument used in the EPOXI mission. Therefore, a realistic earth system model contributed to the direct and quantitative comparison between the observation images and defocus-derived simulation images.

Zhuo and Sim’s defocus mapping tool was used to estimate the defocus [65]. The first row in Fig. 8 shows an observation and six simulation images that were tested with this tool. The six cases of defocus simulation involved varying the distance between the primary mirror and secondary mirror from the diffraction-limited case (0 μm) to 100 μm longer. First, the tool found the location of the edge and made edge maps with the simulation images (the images in the second row in Fig. 8). Then, the tool measured the sharpness of the edges and estimated the blur or defocus maps of the input images. The third row of Fig. 8 represents the defocus maps, in which a high black intensity indicates increasing sharpness or a detailed region, whereas a high white intensity indicates a more blurred region or decreasing sharpness compared to the detailed region. The images in the last row of Fig. 8 show subtraction of the defocus maps of the simulation image from the map of the observation image. As the figure shows, the right edge of the Earth image shifts gradually from black to white as the defocus increases. We estimate that a distance of approximately 60 μm is similar to the defocus of the observation image, whose intensity difference map instantly shifts from a negative to a positive slope.

 

Fig. 8 RGB composite images, edge map, defocus map, and magnified difference map of the observation and simulation maps of defocus aberration.

Download Full Size | PDF

In addition, the slope of the edge spread function (ESF) [66] of this region was measured manually for specific estimates. A 4 × 4 pixel region of the Earth image was selected to calculate the slope of the ESF (as shown by the red box in the fourth row of Fig. 8). Then, the linear slope of the ESF was measured in this pixel region and compared to the slope of the ESF of the observation image. As a result, when the distance between the primary and secondary mirror changed by 55.498 μm from the diffraction-limited condition, the difference in the edge slope was minimized, which indicated the best estimation of the defocused imaging condition for the HRI-VIS instrument. This determination of the defocused imaging condition indicates that 80% of the size of the encircled energy (EE) changes from 22.9 μm to 72.8 μm at full width at half maximum of the PSF from 17.0 μm (0.81 pixels) to 76.8 μm (3.65 pixels) and that the size of the PSF ring changes from 80.0 μm (3.8 pixels) to 190.4 μm (9.1 pixels), based on analysis of the optical design tool comparing the limited performance of diffraction.

Subsequently, we determined that the position of the focal plane should shift to 5.813 mm shorter than its designed position to achieve the best focus. According to previous studies on the HRI-VIS defocus issue, the shift of the focal plane was estimated to be several millimeters shorter than the design position after the launch of the spacecraft [62]. After baking out to reduce the defocus issue, the size of the PSF was estimated to change from 11.6 pixels to 9.0 pixels [62]. This provides good evidence that the quantitative optical performance estimation (9.1 pixels) of the IRT method works correctly, based on direct comparisons between simulated and observed images.

4.4 Disk-averaged spectra

The earth system model produced using the IRT method not only generates a simulated image but also simultaneously provides the line-by-line spectrum of the Earth, unlike other models. The red line in Fig. 9(a) shows the time- and disk-averaged synthetic spectrum for the case of a 24-h simulation from the line-by-line ray tracing computation using the earth system model with the IRT method. The blue line compares the spectra from Robinson et al.’s research with this simulation [24].

 

Fig. 9 Disk-averaged spectra of line-by-line simulation with IRT method.

Download Full Size | PDF

Both graphs indicate the spectral features of the Earth; i.e., its Rayleigh scattering, the absorption lines of oxygen, water vapor, carbon dioxide, and other molecules, and its red edge (0.75–0.95 μm) signature [67]. These spectral features change over time, as shown in Fig. 9(b), which represents a two-dimensional (2-D) spectrum. In Fig. 9(b), the x-axis represents the wavelength, the y-axis represents the observation time, and the color of the 2-D spectrum represents the intensity of the difference in the spectral signature, based on the mean spectrum in Fig. 9(a).

The map on the right of the 2-D spectra corresponds to the viewing geometry of the Earth [68] from the line of sight of the instrument. When the Earth is shown in an ocean-dominant region at times from 13 to 25 h, the albedo of the Earth is darker than in other cases, and the Rayleigh scattering and reflectance of the water’s surface is significant in the blue wavelength region (less than 0.7 μm). On the other hand, a strong red edge signal is shown for the vegetation-dominant region at times from 3 to 12 h in Fig. 9(b).

5. Discussion

This paper presents new techniques for calculating error when comparing observation and simulation data. We were able to detect differences in specific cloud distributions between simulated and observed images. However, a thick cloud distribution, which has a high fraction of clouds, is similar in the case of both the simulation and observation. In this case, it is sufficient to compare simulated and observed images directly. Additionally, this difference does not primarily affect the edge slope estimation that is used for quantitatively evaluating the defocus aberration of the HRI-VIS instrument.

The effective phase function for cloud scattering was derived using the equation presented by Pfeiffer and Chapman for multiple scattering estimation [46]. Piskozub and McKee found that this estimation agrees well with the results of Monte Carlo multiple-scattering calculations [47].

However, the scattering feature, in which the observed image is brighter than the simulation image near specular regions because of the sun glint effect on the Earth’s image, is illustrated in the difference maps in Fig. 8. In future research, the precise backscattering effect produced by the clouds’ surfaces needs to be considered in estimating the spatial radiometric distribution of the Earth.

In this study, detector specifications were only used with quantum efficiency and pixel size, not with radiometric calibration parameters, i.e., pixel response, noise, scattered light, and frame transfer smears [62]. Fortunately, the observation data provided both uncalibrated and calibrated images [37]. The calibrated images were used in the simulation to evaluate the “optical” performance, excluding the electric and radiometric effects of the detector. In addition, the spacecraft used in the Deep Impact/EPOXI mission had a large pointing error that required manual translation of each three-wavelength-band image for creation of the RGB composite. However, there is a chance that residual error remained because of the radiometric calibrations and jitters. In future research, additional calibrations will be studied and applied to estimation of detector performance.

In the calculations performed in this study, the ray tracing algorithm used a geometric optics environment, not a wave optics environment that calculates beam propagation. The diffraction effect was ignored in the imaging simulation because the instrument model had no diffraction-limited condition due to the defocus aberration and because the measurement target was not a pointing source. In addition, the diffraction performance PSF of the evaluated instrument design was confirmed using other optical design tools and matched the results of PSF calibration experiments well.

The ray tracing fluctuation issue has been considered in previous studies [29–31,69]. In this study, the number of rays was greater than 9,000,000 for the imaging simulation and 1,000 for the monochromatic simulation with line-by-line disk-averaged spectra. The simulation with this number of rays kept the radiometric fluctuations below 1% of the calculated image and spectra. For increased calculation efficiency, three cores of a 3.0-GHz quad-core central processing unit (CPU) were used for the imaging simulation, and the other core was used for the disk-averaged spectra simulation for the 20-h per simulation case.

This paper proposes a new global earth system and an instrument model for use in a ray tracing environment. The quantitative and reasonable optical performance of the Earth-observing instrument with the IRT method was verified using HRI-VIS the Earth observation data and spectral results from other studies. This earth system model, which includes instrument consideration, is suitable for use as an end-to-end simulation tool that is interchangeable and combines models in a unified environment. The earth system model takes into account most of the sub-models used to constitute the Earth. This model was constructed using time-, spatial-, and spectral-resolved input data to obtain the time- and disk-averaged signature of the Earth. The instrument model was able to produce simultaneous imaging and spectroscopic results from line-by-line ray tracing. The optical performance estimation was expanded using optical design variation in the ray tracing method.

6. Concluding remarks

In this study, a new global earth system model is introduced for the optical performance evaluation of various space instruments. The model is unique in that it provides simultaneous imaging and spectroscopic results, which offers high compatibility with various imaging, spectroscopic, and hyperspectral imaging instruments.

The global earth system model also provides fully resolved spatial, spectral, and temporal coverage of the Earth and includes sub-models for the molecular atmosphere, clouds, aerosols, land surfaces, snow cover, sea ice, and ocean surfaces. The cloud sub-model alone is capable of evaluating performance when observing a cloudy sky.

The simulation results for the image- and disk-averaged spectra obtained with the proposed earth system model were compared with observational data from HRI-VIS measurements from the EPOXI mission. The defocus aberration of HRI-VIS was verified using performance estimation via direct comparison of simulated and observed images. The results of this study shows that the combination of the proposed global earth system model and IRT is a powerful tool to simultaneously evaluate the imaging, radiometric, and spectral performance of instruments used in observing the Earth. End-to-end simulations performed using this model will lead to efficient progress toward the achievement of the scientific goals of optical instrument missions.