Atmosphere and Terrain Coupling Simulation Framework for High-Resolution Visible-Thermal Spectral Imaging over Heterogeneous Land Surface

: Realistic modeling of high-resolution earth radiation signals in the visible-thermal spectral domain remains difﬁcult, due to the complex radiation interdependence induced by the heterogeneous and rugged features of land surface. To ﬁnd the trade-off between accuracy and efﬁciency for image simulation, this paper established a uniﬁed simulation framework for the entire visible-thermal spectral domain, based on the energy balance between solar-reﬂected and thermal radiation components over rugged surfaces. Considering the joint contributions of atmospheric and topographic adjacency effects, three spatial–spectral convolution kernels were uniformly designed to quantify the topographic irradiance, the trapping effect, and the atmospheric adjacency effect. Radiation signal simulation was implemented in three forms: land surface temperature (LST), bottom of atmosphere (BOA) radiance, and top of atmosphere (TOA) radiance. The accuracy was validated with onboard data from China’s Gaofen-5 visual and infrared multispectral sensor (VIMS) over rugged desert. The simulation results demonstrate that the root mean square of relative deviations between the simulated and onboard TOA radiance are related to terrain, as 3–17% and 6–38% for the summer and winter scene, respectively. The evaluation of radiance components indicates the utility of the simulation framework to quantify the uncertainty associated with atmosphere and terrain coupling effects, in the sensor design and operation stages. adjacency irradiance contribution to L , and introducing extra spatial varieties to the upward transfer path from different points of the ground scene to the sensor. As shown in Figure 1, the spectral signal recorded by an imaging sensor is the sum of seven interaction components derived from solar and thermal radiations. topographic irradiance, trapping effect, and atmospheric adjacency effect with three spatial-spectral convolution kernels related to both atmospheric and terrain features, avoiding the commonly used assumption of independence between atmospheric and topographic adjacency effects; (II) it provides a novel method to establish layer atmosphere coefﬁcients look-up tables based on multiple output ﬁles of MODTRAN, which beneﬁts the imaging simulation for various viewing altitude and angles. The simulation framework, and the produced data, could be used for the performance predictions of earth-observing systems, and serve to optimize spectral imaging missions. The accuracy of the TOA radiance simulation was validated using China’s Gaofen-5 VIMS data. The simulation experiment indicates that the proposed simulation framework is reliable. The spectral root mean square of relative deviations between simulated radiance and onboard radiance is less (larger) than 17% (3%) and 38% (6%) in the summer and winter scene, respectively. Some statistical evaluations of the atmospheric and topographic effects were also implemented. It indicates that the uncertainty associated with atmosphere and terrain coupling effects should be quantiﬁed under the whole cycle, from the concept design stage to the application stage, for a spectral imaging mission of earth.


Introduction
High-resolution spectral imaging systems in the visible to thermal spectral domain play an essential role in the quantitative remote sensing of land surface [1,2]. Generally, the visible and short-wave infrared (VSWIR) and the thermal infrared (TIR) imaging spectrometers are independently operated to gather remote sensing datasets. However, the simultaneous measurement of longwave and shortwave features ensure consistency across both ranges, and provide vantages on time-varying application scenarios, such as ecosystem investigations and urban mapping. The recent study of Fahlen et al. [3] proves that a consistent retrieval of visible to thermal data reduces the posterior uncertainties for the surface and atmosphere features, such as the land surface temperature (LST) and air temperatures. In recent years, China launched the Gaofen-5 series of satellites [4][5][6] to gather hyperspectral images in the 0.4-2.5 µm and multispectral images in the 0.45-12.5 µm spectral domain. China started the design and research of the geostationary orbiting hyperspectral imager over the visible to thermal spectral domain. Continuous spectral imaging over the entire visible-thermal domain, with a ground sampling distance (GSD) scattering approximation were adopted in the quantification processes of atmospheric and terrain adjacency effects, resulting in three spatial-spectral convolution kernels. Molecular/aerosol scattering phase functions, topographic features, surface heat balance, and ground heterogeneity were taken into consideration in the simulation framework. A particular atmosphere parameter look-up table based on MODTRAN 5 was also introduced. Section 3 is devoted to the experimental results. The simulated radiance images were validated using China's Gaofen-5 visual and infrared multispectral sensor (VIMS) [6] dataset. The atmospheric and topographic contributions on radiation signal were then evaluated and compared with the VIMS' sensor performance. Lastly, the discussions and conclusions are provided in Sections 4 and 5, respectively.

Signal Component of At-Sensor Radiance
The measurement function of a radiance signal within a sensor's linear responding region was properly simplified to a convolution process, with the abstract conception of sensor equivalent response function R sensor , and a random noise process with the equivalent radiance noise L n , as shown in (1). Imaging the atmosphere as a medium layer between the surface and sensor, the abstract conception of atmospheric equivalent response function R atmo , also described the atmospheric direct and diffuse influence on radiance signals.
L res_TOA = L TOA ⊗ R sensor + L n = (L GND ⊗ R atmo + L a ) ⊗ R sensor + L n (1) where L res_TOA and L TOA are the restored onboard and the truth TOA radiance; L GND is the ground-leaving radiance carrying features of surface reflectivity, emissivity, and temperature; and L a is the atmosphere path radiance without ground participations. For flat land surfaces, the atmospheric equivalent response function mainly described the integration of the atmospheric attenuation of ground signals and the atmospheric adjacency effect. For the rugged land surface investigated in this paper, topography also affected R atmo , by impacting solar and sky downward irradiance on the surface, increasing the adjacency irradiance contribution to L GND , and introducing extra spatial varieties to the upward transfer path from different points of the ground scene to the sensor. As shown in Figure 1, the spectral signal recorded by an imaging sensor is the sum of seven interaction components derived from solar and thermal radiations. and terrain adjacency effects, resulting in three spatial-spectral convolution ker Molecular/aerosol scattering phase functions, topographic features, surface heat bala and ground heterogeneity were taken into consideration in the simulation framewor particular atmosphere parameter look-up table based on MODTRAN 5 was introduced. Section 3 is devoted to the experimental results. The simulated radi images were validated using China's Gaofen-5 visual and infrared multispectral se (VIMS) [6] dataset. The atmospheric and topographic contributions on radiation si were then evaluated and compared with the VIMS' sensor performance. Lastly discussions and conclusions are provided in Section 4 and Section 5, respectively.

Signal Component of At-Sensor Radiance
The measurement function of a radiance signal within a sensor's linear respon region was properly simplified to a convolution process, with the abstract conceptio sensor equivalent response function sensor R , and a random noise process with equivalent radiance noise n L , as shown in (1). Imaging the atmosphere as a medium l between the surface and sensor, the abstract conception of atmospheric equiva response function atmo R , also described the atmospheric direct and diffuse influenc radiance signals. For flat land surfaces, the atmospheric equivalent response function ma described the integration of the atmospheric attenuation of ground signals and atmospheric adjacency effect. For the rugged land surface investigated in this pa topography also affected atmo R , by impacting solar and sky downward irradiance on surface, increasing the adjacency irradiance contribution to GND L , and introducing e spatial varieties to the upward transfer path from different points of the ground scen the sensor. As shown in Figure 1, the spectral signal recorded by an imaging sensor i sum of seven interaction components derived from solar and thermal radiations. Component 1: spectral radiation originating from the thermal emission of the rug land surface, after upward propagating in the atmosphere, finally reaches the sen Figure 1. Schematic of the sun-surface-atmosphere-sensor system radiation signal composition. Seven components of at-sensor radiance are numbered sequentially, and shown as 1 to 7. Component 1: spectral radiation originating from the thermal emission of the rugged land surface, after upward propagating in the atmosphere, finally reaches the sensor. Component 2: spectral radiation originating from the solar direct irradiance, after reflecting by land surface and upward propagating in the atmosphere, finally reaches the sensor. Component 3: spectral radiation originating from the scattered solar radiation and the thermal emission of atmosphere, after reflecting by land surface and upward propagating in the atmosphere, finally reaches the sensor. Component 4: spectral radiation originating Remote Sens. 2022, 14, 2043 4 of 30 from the adjacent irradiation of the rugged land surface, after reflecting back from the surface and upward propagating in the atmosphere, finally reaches the sensor. Component 5: spectral radiation leaving from the adjacent pixels, containing all the contributions of component 1 to 4, is successively trapped by the atmosphere, reflected by land surface, then propagated in the atmosphere, and finally reaches the sensor. Component 6: spectral radiation originating from the exitance of adjacent pixels, containing the contributions of surface thermal emissions and reflections, after upward scatting by atmosphere, finally reaches the sensor. Component 7: spectral radiation originating from the atmospheric scattering of solar energy and the thermal emission of atmosphere, is partly scattered into the IFOV of the sensor, without reaching the land surface. Among these components, components 1 to 6 introduced the influence of land surface to the total TOA signal. Components 1 to 3 describe the direct emissions and reflections of target pixel P, while components 4 to 6 describe the radiation exchanges among target, background, and atmosphere.
As expressed in (2), calculations of the seven components were modularized and summed to the TOA radiance signal, where three convolution processes were introduced to describe the adjacency effects. The convolution kernel F was used to describe the spatial and spectral distribution of the background radiation contribution ratio caused by atmospheric adjacency effect. The convolution kernel G was used to describe the spatial and spectral distribution ratio of adjacent terrain irradiance reaching the target P. The convolution kernel H was used to calculate the background equivalent reflectance that affects the successive reflection irradiance between ground and atmosphere. Compared to our previous work on solar-reflected spectral domain [36], the thermal emission radiance of land surface (εL B ) was introduced as component 1, and the thermal diffused irradiance from sky (E in_therm ) was introduced in component 3. The two sources of thermal radiation then provided joint increase to components 4, 5, and 6.
where L TOA_rug and L GND_rug are TOA and BOA radiance over the rugged land surface, respectively; τ dir and τ di f are the direct and diffuse transmittance matrix, respectively; ε is the surface emissivity; ρ is the surface reflectance; L B is the blackbody radiance, which is calculated using the Plank equation and LST; E in_dir , E in_di f , and E in_therm are the direct solar irradiance on flat ground, the diffuse solar irradiance on flat ground, and the diffused thermal irradiance on the flat ground, respectively; SF is the shade factor matrix with each element spread from 0 to 1, describing the cast shadow; θ i is the solar illumination angle (SIA, angle between the surface normal and the sun ray), describing the self-shadow effect; V sky is defined as the relative proportion of solid angle towards the sky; S is the hemispherical scattering albedo of BOA; and E 1 out_rug is the equivalent initial radiation of the trapping effect.
To simulate the components related to surface emission over rugged land surface, the modeling of heterogeneous high-resolution LST was necessary. Instead of using the analysis-ready LST data from satellite observations, as adopted in previous works [33], we tentatively combined the above simulation of the slope-leaving radiation process with the surface heat balance theory, to simulate the LST over rugged land surface. This strategy provided sub-pixel LST scene with a GSD smaller than the analysis-ready LST data. Another advantage of this strategy was that it avoided the correlativity of other sensors being additionally introduced to the sensor performance prediction tasks. As auxiliary products of the simulation framework, the surface net radiation flux and the flux caused by topographic irradiance were generated by integrating the corresponding spectral flux in the spectral range of 0.2 µm to 50 µm. The surface heat balance function (10) was solved using the newton iteration method to simulate LST of the rugged scene.
where R n is the surface net radiation flux, defined as the pixel radiation net flux caused by E in_rug − E 0 out_rug ; F adj is the net flux caused by the topographic irradiance; HS is the sensible heat flux; LE is the latent heat flux; and HC is the heat conduction flux.

Atmospheric and Terrain Effects Modeling
The topographic irradiance only exists due to the relief around the target. Figure 2 illustrates the geometrical relationship between the target and the adjacent pixel that illuminates the target. The rugged surface was assumed as a diffuse irradiance source. The convolution kernel G describes the ratio of the radiation of each surrounding pixel M that reached the target pixel P. According to the solid angle projection law, the receiving irradiance of facet P could be expressed as (10).
where L M is the equal radiative of the lambert radiation source M; θ P is the angle between the line MP and the normal of P; and E out, M is the slope-leaving radiation of facet M, calculated using (4); dω M→P = cos θ M dA M /l MP 2 is solid angle of facet M towards P, A M is the area of M, and l MP is the distance between P and M. the flux caused by topographic irradiance were generated by integrat corresponding spectral flux in the spectral range of 0.2 μm to 50 μm. The surf balance function (10) was solved using the newton iteration method to simulat the rugged scene. The topographic irradiance only exists due to the relief around the target. illustrates the geometrical relationship between the target and the adjacent p illuminates the target. The rugged surface was assumed as a diffuse irradiance sou convolution kernel G describes the ratio of the radiation of each surrounding that reached the target pixel P. According to the solid angle projection law, the r irradiance of facet P could be expressed as (10).  Considering the spectral transmittance of the path from M to P, the visibility P, and the slope angle of M, the convolution function G is expressed in (11). G only a function of distance, but also wavelength, because the path transmittanc with wavelength. As the optical distance between M and P increased, the top irradiance decreased. Thus, the spatial range of G was determined by the terrain transmittance. The convolution function was sampled along the orthogonal co system of the raster scene data.  Considering the spectral transmittance of the path from M to P, the visibility of M by P, and the slope angle of M, the convolution function G is expressed in (11). G was not only a function of distance, but also wavelength, because the path transmittance varied with wavelength. As the optical distance between M and P increased, the topographic irradiance decreased. Thus, the spatial range of G was determined by the terrain and path transmittance. The convolution function was sampled along the orthogonal coordinate system of the raster scene data.
where k MP is the equivalent extinction coefficient along the M to P path; V M→P is the visibility factor describing whether the irradiance could reach P along the M to P path, where 1 represents visible and 0 represents invisible; and α M is the slope angle of pixel M. Pixel P was also the diffuse irradiance source to M. The energy balance of topographic irradiance was solved using iteration algorithms [37]. Generally, three or four times of iteration made the energy balance function converge, even for deep-shaded areas. As a trade-off between accuracy and calculation cost, the global convolution calculation was limited to three times in this study. Figure 3 illustrates that the radiation source of the trapping irradiance is the spatial average of the ground-leaving radiation, which included the target emission, as well as the target reflection of the downward irradiance and topographic irradiances. Only the upward potions of ground-leaving radiation within the solid angle (Ω sky ) toward sky formed the trapping irradiance. Over a rugged scene, irradiance at ground varied from pixel to pixel. The coupling irradiance caused by trapping effect could also be approached using the Monte Carlo method [38,39], but it is an extremely time-consuming process for spectral image simulation. Over a flat scene, the successive reflections and scattering between the surface and the atmosphere, the so-called trapping effect, was generally modeled using the multiply operation between the ground-leaving radiance, with a gain term of 1/(1 − ρ e S) [40]. With the increase in scattering times, the contribution of the trapping effect decreased gradually. We used the gain term concept, considered the terrain influence on the adjacency average reflectance (ρ e_rug ), and the average trapping irradiance (E out_rug ) over a rugged scene. The trapping factor over a rugged scene was defined as (12). Trap f actor(x, y, λ)

Trapping Effect
Assuming the ground target had Lambertian reflectance characteristics, the upward portion ratio of ground-leaving radiation was also described by the sky view factor. The spatial average of upward radiation was calculated as (13).
where the weighted integral of slope-leaving radiation is conducted within an adjacent range around each target pixel. As the aerosol scattering was the main origin of atmospheric hemispheric albedo, the integral calculation was operated within the same spatial range with atmospheric adjacency effect.
Taking the successive reflection process as the arithmetic of geometric sequence, the common ratio was calculated using the ratio between the two sequential upward radiation from the rugged surface, as shown in (14). Therefore, the trapping convolution kernel function was described as (15).

Atmospheric Adjacency Effect
The surface-leaving radiation escaped the shackles of land surface and was partially transmitted by the atmosphere into the sensor's IFOV, and partially scattered into the adjacent IFOVs. The scattered portion, the so-called atmospheric adjacency effect, was generally modeled using the convolution of ground-leaving radiation and an atmospheric PSF [41]. The shape of the atmospheric PSF depended on the aerosol scattering phase function and aerosol optical depth [32]. For a rugged surface, the topography actually changed the atmospheric scattering actions. In our simulation flow, the single scattering method [12,36,42] was introduced to describe the spatial contribution ratio of upward leaving-radiation from adjacency pixel to the atmospheric adjacency effect. The atmospheric diffuse transmission was introduced to describe the influence of the multiple scattering of aerosols and molecules. Figure 4 is a schematic diagram of single scattering geometry. The atmosphere was vertically stratified to N layers, according to altitude, and it was assumed that the absorption and scattering characteristics of each layer were homogeneous. The photon emitted by an adjacency pixel M was scattered once by a layer to the viewing direction. The sum of the scattered portions by the N layers described the atmospheric PSF (PSF atmo ) in (16).
where ς i is the scattering coefficient of the ith layer; γ i is the scattering angle, i.e., the angle between the path from the M to the atmospheric scatter point Q i and the viewing direction; P(γ i ) is the scattering phase function; ω Q i M is the angle between surface normal direction and the path from M to Q i ; τ MQ i and τ Q i Z v are transmittances of the path from M to Q i and the path from Q i to the sensor, respectively, calculated using the Beer-Bouguer-Lambert law [43]; l MQ i is the distance between M and Q i ; and θ v is the sensor viewing angle for the slope pixel. The distance and intersection angles between lines in Figure 4 were determined by geometrical relations among ground pixels, atmospheric layers, and the sensor. The scattering coefficient was a vertical varying atmospheric parameter, while the scattering phase function and path transmittance were functions of both atmosphere and terrain parameters. Thus, the atmospheric PSF varied with atmospheric condition, pixel location, imaging time, viewing direction, adjacency distance, and wavelength. Both the scattering of molecules and aerosol particles were taken into consideration in our simulation framework. The Rayleigh scattering phase function was used to describe molecule scattering intensity distribution along with angle, as shown in (17). For the simplification of Mie scattering by aerosol, the Henyey-Greenstein phase function [43] was introduced, as shown in (18), with an asymmetry factor (g) varying from −1 (completely back scattering) to 1 (completely fore scattering).
Remote Sens. 2022, 14, 2043 10 of 30 The relative weight of each adjacency pixel was calculated by normalizing the atmospheric PSF, as shown in (19). The relative weight function was then used as the convolution kernel function in our simulation framework.

Atmospheric and Topographic Parameter Generation
The above-mentioned calculations of atmosphere and terrain coupling effects were all related to topographic parameters. The slope, aspect, and sky view factor of natural bare earth surface are intrinsic characteristics of topography, and do not change considerably over time. We used the digital elevation model (DEM) to provide the basic data for descriptions of the surface relief features. The topographic parameters were gathered using the method shown in Table 1. In order to ensure the accuracy of the simulation, the spatial resolution of DEM was higher than the simulated imaging sensor's resolution.
The nine atmospheric coefficients in Table 2 are essential parameters for the calculations of TOA and BOA radiance, and were derived from MODTRAN output files. The MODTRAN interrogation technique used in other works invoked the MODTRAN under different surface albedos, and recalculated approximate coefficients of the whole atmospheric medium [31][32][33][34]. The zero surface albedo conditions of these methods leads to miscalculations when extracting parameters from MODTRAN [44]. Different from these previous works, we avoided the recalculations with different surface albedos and gathered layer coefficients of atmosphere. The "fort69", "flx", "acd", and "tp7" files, which are exported by MODTRAN 5.3.2 [16], were applied for atmospheric coefficient generation in this paper. Calculation methods are also described in Table 2. The option of "IEMSCT" in Card1 of MODTRAN gave the separate simulation of solar diffuse and thermal diffuse irradiance. Considering the variations of atmospheric coefficients due to surface elevation over imaging scene, MODTRAN was invocated at different surface elevation conditions for coefficients interpolation. The spectral interval option of MODTRAN was set 1 cm −1 in our simulation framework. The atmospheric parameters were convoluted using the sensor's spectral response function (SRF). Table 1. Terrain parameters generation.

Terrain Parameter Symbol Calculation Method
Slope angle α α = arctan f x 2 + f y 2 where the inverse square distance weighted spatial interpolation of DEM [45], was used to calculate, f x and f y along a west-east and south-north direction, respectively.

Aspect angle
Shade factor SF Shade factor was calculated using the depth buffer method.

Sky view Factor
where h φ is defined as the angle between the zenith and the horizon surface at a given orientation direction φ [46,47]. Table 2. Atmospheric parameters generation.

Atmospheric Coefficients Symbol Calculation Using the MODTRAN Output Data
Hemispherical scattering albedo of BOA 1 S S = interp(DEM, SPHERICAL ALBEDO FROM GND) where z layer,i is the height interval of the ith layer.

Scattering coefficient of layers
The coefficient is derived from information in the "acd" file. 2 The coefficient is derived from information in the "tp7" file. 3 The coefficient is calculated using information in the "flx" file. 4 The coefficient is calculated using information in the "fort69" file.

Materials for Simulation Case Study
As mentioned in Sections 2.1 and 2.2, we emphatically introduced a simulation methodology of the atmosphere and terrain coupling effects over rugged surfaces. The DEM and reflectance were general inputs of the simulation framework, as well as the two most variable terms to describe the land surface conditions. The final output of the simulation framework was the TOA radiance signal. In this paper, a simulation experiment was designed to validate the simulation fidelity of the time-varying TOA radiance over rugged surfaces, by comparing the simulated and the onboard radiance data.
The spaceborne data gathered by the VIMS of China's Gaofen-5 satellite were selected as reference for simulation validation. The Gaofen-5 satellite, the fifth satellite of China's high-resolution Earth observation system (CHEOS), was launched on May 9, 2018, to an orbit of 705 km altitude. There are two high-resolution spectral imaging sensors onboard the Gaofen-5 satellite. The advanced hyperspectral imager (AHSI) is observing the Earth with 330 spectral bands in the visual to short-wave infrared spectral domain (0.4-2.5 µm), with a GSD of 30 m. The VIMS provide 12 separate bands of images in the visual to thermal domain. A 20 by 15 km area of Mingsha mountain in Dunhuang, Gansu Province of China, was selected as the study area. This is a rugged sandy mountain, with an altitude ranging from 1100 m to 1700m. The central point coordinates are 40 • 1 34.78"N and 94 • 37 19.96"E. Two imaging data, gathered on 14 July 2019 and 12 December 2019, were selected to represent two typical conditions of the different environments of summer and winter, respectively. A brief introduction of the imaging scenes and the instrument performance information of VIMS are shown in Table 3. Table 3. Imaging information and instrument performance parameters of Gaofen-5 VIMS.
The simulation and evaluation flow are shown in Figure 5. The overall process was implemented by the following steps: (1) the calculation of topographic parameters using DEM data; (2) the calculation of atmospheric parameter look-up tables using MODTRAN; (3) the generation of the hyperspectral reflectance scene in the visual to thermal domain, using a linear mixing model of reflectance spectra and the abundance map derived from AHSI data; (4) the generation of the LST scene based on the heat balance model of rugged surfaces; (5) the calculation of the rugged irradiance from sun and sky, the topographic irradiance effect, the trapping effects, and the atmospheric adjacency effects in sequence to generate the TOA radiance; and (6) the radiometric calibration of the VIMS onboard data, and the comparison between the simulated and restored TOA radiance images.
The simulation and evaluation flow are shown in Figure 5. The overall process was implemented by the following steps: (1) the calculation of topographic parameters using DEM data; (2) the calculation of atmospheric parameter look-up tables using MODTRAN; (3) the generation of the hyperspectral reflectance scene in the visual to thermal domain, using a linear mixing model of reflectance spectra and the abundance map derived from AHSI data; (4) the generation of the LST scene based on the heat balance model of rugged surfaces; (5) the calculation of the rugged irradiance from sun and sky, the topographic irradiance effect, the trapping effects, and the atmospheric adjacency effects in sequence to generate the TOA radiance; and (6) the radiometric calibration of the VIMS onboard data, and the comparison between the simulated and restored TOA radiance images.
The simulation and evaluation flow are shown in Figure 5. The overall process was implemented by the following steps: (1) the calculation of topographic parameters using DEM data; (2) the calculation of atmospheric parameter look-up tables using MODTRAN; (3) the generation of the hyperspectral reflectance scene in the visual to thermal domain, using a linear mixing model of reflectance spectra and the abundance map derived from AHSI data; (4) the generation of the LST scene based on the heat balance model of rugged surfaces; (5) the calculation of the rugged irradiance from sun and sky, the topographic irradiance effect, the trapping effects, and the atmospheric adjacency effects in sequence to generate the TOA radiance; and (6) the radiometric calibration of the VIMS onboard data, and the comparison between the simulated and restored TOA radiance images.
The simulation and evaluation flow are shown in Figure 5. The overall process was implemented by the following steps: (1) the calculation of topographic parameters using DEM data; (2) the calculation of atmospheric parameter look-up tables using MODTRAN; (3) the generation of the hyperspectral reflectance scene in the visual to thermal domain, using a linear mixing model of reflectance spectra and the abundance map derived from AHSI data; (4) the generation of the LST scene based on the heat balance model of rugged surfaces; (5) the calculation of the rugged irradiance from sun and sky, the topographic irradiance effect, the trapping effects, and the atmospheric adjacency effects in sequence to generate the TOA radiance; and (6) the radiometric calibration of the VIMS onboard data, and the comparison between the simulated and restored TOA radiance images.    The necessary auxiliary data of the simulation experiment are shown in Figure 6. The land cover map of the study area was taken from the sentinel 2 dataset of the European Space Agency [48], with a GSD of 10 m, as shown in Figure 6a. It illustrates that there are mainly four classes of land cover (including cropland, build-up, desert, and water) in the specified study area. The flat vegetation and urban areas are located in the northwest (upper right) corner of the study area. The Mingsha mountain area is a rugged and uniform sandy mountain, occupying the other areas of the scene. The ALOS DEM dataset [49], with a GSD of 12.5 m, as shown in Figure 6b, was used to extract the subpixel topographic features. The GSD of the VIMS and AHSI data are all about 30 m. Thus, the pixels of the flat vegetation and urban areas are generally a mixture of several endmembers. The reflectance product of the synchronous AHSI data gave the reflectance cube of the mixed pixels in the 0.4-2.5 µm, while there are no suitable hyperspectral imaging data to produce the reflectance of the mixed pixels.
The linear mixing model was applied to establish the reflectance of each mixed pixel. The AHSI reflectance products were resampled to the scene of 10 m, and registered to the land cover map. The retrieved reflectance spectra of the four typical land covers were extracted from the region of interests (ROI) of the resampled AHSI reflectance products. The ROI of the desert (sand) was marked in the flat side of the Mingsha mountain, near the urban area. As the rugged Mingsha mountain is a uniform sandy mountain, we assigned these pixels with the same reflectance spectra of the desert (sand) in the artificial scene. The mean spectra of these retrieved reflectance were taken as the four typical endmembers, and are shown in Figure 6c. An abundance map, with a GSD of 10 m, was gathered using the AHSI data and the endmember spectra. The abundance map was then used to establish the infrared spectral domain, by mixing the infrared reflectance spectra of the four kinds of endmembers. The infrared endmember spectra were gathered from field measurement results and the spectral library. The ground-leaving spectra of sand on Mingsha mountain, river water, and the flat Gobi Desert were measured using a D&P 102F FTIR spectrometer in July of 2018, then inverted to emissivity using the previously published TES method [50]. We converted these retrieved land surface emissivity (LSE) to reflectance spectra using the Kirchhoff's law [51] for the infrared scene simulation, and the vicarious calibration of onboard data. The reflectance spectra of the cropland, built-up of 3 to 14 µm, were replaced using reflectance spectra from the ASTER spectral library [52]. The mean spectrum of soil and vegetation was taken as the reflectance of the cropland endmember. A spectrum of construction in the spectral library was taken as the build-up endmember.
Another variable factor that had an essential influence on the radiance signals was the local atmospheric conditions when gathering the imaging data. The radar observation of relative humidity profiles, shown in Figure 6e, were used to represent the water vapor conditions of the troposphere on the two experimental days. Temperature profiles generated from the fifth generation ECMWF atmospheric reanalysis of the global climate (ERA5) dataset [53], and in-site measurements shown in Figure 6f, were averaged to represent the atmosphere temperature under the elevation of 40 km. These profiles were resampled and jointed with the standard profile of "1976 US Standard Atmosphere" mode in MODTRAN for higher layers, to describe the condition of the whole atmosphere medium under 100 km.
As the validation of the performance modeling of instruments was not the major task of this paper, a simplified sensor degeneration simulation of VIMS was carried out. All the atmospheric coefficients and spectra were resampled to spectral variates, using the SRF shown in Figure 6c,d. The SRFs of B1 to B7 are described using the spectral band range provided in the VIMS L1 dataset shown in Table 1, based on Gaussian functions, while the SRFs of B8 to B12 used here were the SRF provided by the numerical weather prediction satellite application facility [54].  The designed values of instrument performance parameters including GSD, the signal to noise ratio (SNR), the noise equivalent temperature difference (NE∆T), and the modulation transfer function (MTF), shown in Table 3, were used to describe the sensor degeneration effects on the TOA signal. To simulate the sensor's spatial performance, the MTF at Nyquist sampling frequency was used to establish a spatial convolution function, as show in (20) and (21). In order to avoid the direction error and limit the extra aliasing effect in simulated data, the spatial interpolation-resampling technology [55] was applied to the simulated TOA radiance. The simulated images were spatially resampled along and across orbit direction, using the rational polynomial coefficients (RPC) provided by the VIMS L1 dataset. Each imaging pixel was interpolated into sub-pixels, then convoluted with the spatial filter function in (20), and the radiance of sub-pixels were merged to the imaging pixel area of VIMS. The noise equivalent radiance (NER) was generally used to describe the standard deviation of the sensor random noises on TOA radiance [7]. We used the designed SNR and NE∆T to calculate the NER of each spectral band of VIMS, as shown in (22). The random noise of VIMS was simplified as white Gaussian noise with NER as standard deviation, and added to the TOA images.
where n b is band number of VIMS; T n is the blackbody temperature that NE∆T is defined; and ∆T is a temperature difference for the linearization of Plank function.
Considering the unavailability of the radiometric calibration coefficients of VIMS, the cross calibration or vicarious calibration process is generally recommended [5]. For the two scenes in this paper, the simultaneously acquired hyperspectral data by the welled calibrated hyperspectral imager of AHSI was used in the cross calibration for the B1 to B6 of VIMS. Both the VIMS and AHSI data were orthographically corrected, based on RPC model. The VIMS data was then spatially resampled to a GSD of 30 m. The flat Gobi Desert, near Dunhuang airport, with an area of 9 square km, was selected for the cross calibration between VIMS and AHSI, as shown in Figure 6a. To eliminate the spectral radiation difference caused by the unreliable calibration in the mid-wave to long-wave infrared bands (B7 to B12 of VIMS), two reference points of water and flat Gobi Desert were selected for vicarious calibration, as shown in Figure 6a. The band of B9 was selected to retrieve the LST of the reference points. The reflectance and LST of the two reference points were converted to TOA radiance using the MODTRAN code, and used for the vicarious calibration of the study areas. The restored TOA radiance images after the linear cross calibration and vicarious calibration were used to evaluate the simulated images. The relative deviation (RD) of every pixel, as shown in (23), and the root mean square of relative deviation (RMSRD) for each spectral band, shown in (24), were calculated as the evaluation indexes.
where L SI MU is the simulated TOA radiance; L ONBOARD is the restored TOA radiance from onboard data; and n is the total pixel numbers of rugged sand mountain over the study area.

Validation of the Simulation Results
To the best of our knowledge, the simulated imaging data based on the flat surface assumption were traditionally used as testing dataset for measurement performance prediction and retrieval method evaluation. Under the flat surface assumption, the variable surface characters are only caused by the land cover. The topographic adjacency effects are not considered in these data. Thus, we set the DEM of each pixel as 1379 m, which is the regional average value of the DEM of the study area, to simulate the data under the flat surface assumption. In this section, the simulated data considering both the atmospheric and terrain effects is not only compared with the onboard data, but also with the simulated data under the flat surface assumption.
The simulated LST scenes of the two imaging days are the intermediate results of the proposed simulation framework and shown in Figure 7a-d. The spatial-temporal variant phenomenon of LST is modeled in the proposed simulation framework. The LSTs for the () summer scene on 14 July 2019 are mainly distributed in the range of 280 K to 345 K, while for the winter scene on 21 December 2019, the LSTs are in the range of 260 K to 310 K. The standard deviation of LST under the flat surface assumption is 6.5 K and 1.2 K for the summer and winter scene, respectively, while the standard deviation of the LST under the rugged surface assumption is 7.0 K and 7.3 K for the summer and winter scene, respectively. It is obvious that the flat surface assumption provides less information in the LST signal than the rugged surface assumption. Compared with the product matrixes of the shade factor, and the cosine of SIA shown in Figure 7e,f, it is realistic that the simulated LSTs of the mountain area under the rugged surface assumption increase with solar illumination. Thus, it is reasonable to use the radiative transfer model over rugged surfaces, combined with heat balance theory, to construct of hypothetical LST scenarios. effects are not considered in these data. Thus, we set the DEM of each pixel as 1379 m, which is the regional average value of the DEM of the study area, to simulate the data under the flat surface assumption. In this section, the simulated data considering both the atmospheric and terrain effects is not only compared with the onboard data, but also with the simulated data under the flat surface assumption. The simulated LST scenes of the two imaging days are the intermediate results of the proposed simulation framework and shown in Figure 7a-d. The spatial-temporal variant phenomenon of LST is modeled in the proposed simulation framework. The LSTs for the () summer scene on 14 July 2019 are mainly distributed in the range of 280 K to 345 K, while for the winter scene on 21 December 2019, the LSTs are in the range of 260 K to 310 K. The standard deviation of LST under the flat surface assumption is 6.5 K and 1.2 K for the summer and winter scene, respectively, while the standard deviation of the LST under the rugged surface assumption is 7.0 K and 7.3 K for the summer and winter scene, respectively. It is obvious that the flat surface assumption provides less information in the LST signal than the rugged surface assumption. Compared with the product matrixes of the shade factor, and the cosine of SIA shown in Figure 7e,f, it is realistic that the simulated LSTs of the mountain area under the rugged surface assumption increase with solar illumination. Thus, it is reasonable to use the radiative transfer model over rugged surfaces, combined with heat balance theory, to construct of hypothetical LST scenarios.  The true color image comparisons between the simulated TOA radiance under both the flat and the rugged surface assumption are shown in Figure 8. Figure 8a,b present the true color images of the simulated TOA radiance in summer, and Figure 8c,d present the true color images of the simulated TOA radiance in winter. For the vegetation and urban areas in the winter scenes, the result under the flat surface assumption (Figure 8c) is more realistic than the result under the rugged surface assumption (Figure 8d). The vegetation and urban areas have small slope angles, less than 10 • . Thus, the uncertainty of the topographic parameters increases the error of the topographic irradiance that finally transferred to the TOA radiance. As the solar irradiance is lower in winter than in summer, the TOA uncertainty caused by topographic uncertainties of the flat surface becomes obvious in winter. However, the TOA radiance of the rugged Mingsha mountain area is uniform under the flat surface assumption, because of the absence of topographic effects. The color and shade effects of the simulated TOA radiance under the rugged surface assumption as shown in Figure 8b,d are more consistent with the quick views of onboard images, especially over the rugged areas in Table 3. As the test data for long-term performance prediction, the spatial relation of the light and shade of pixels with different land covers in the simulation image conforms to the actual situations. Visual comparisons of the simulated TOA radiance images under the rugged surface assumption, and the onboard radiance images of all 12 spectral bands are shown in Figure 9  The TOA spectra of several typical areas and pixels, as shown in Figure 8, are compared with the onboard spectra gathered by the Gaofen-5 satellite. Figure 10 shows the TOA spectra comparisons among representative areas with different land covers. The ROI potions of cropland and water are marked in Figure 8b,d. The whole area of the rugged Mingsha mountain is taken as the ROI for the land cover of desert (sand). The mean spectra, with their adjacent interval within the one standard deviation range, are compared and shown in Figure 10a,c,d. One pixel of the road in the urban area is also selected as an example, representing the land cover of construction. In all four of the subfigures, the onboard and the simulated spectra have similar trends in spectral domain. The spectral differences between different imaging times are remarkable for all the four representative areas. In other words, the summer spectra and the winter spectra are distinguished within the one time of standard deviation areas. In addition, the cropland and the urban pixels have more vegetation covers in summer compared to winter. Accordingly, the simulated summer spectra in Figure 10a,b display obvious absorption features in the red spectral band, which is consistent with the onboard data. It is also obvious that the absolute value of the spectra standard deviation in the visual to SWIR bands are larger than the values in the MWIR and LWIR bands. The relative values of standard deviation in the spectral range of 1-5 µm are larger than other bands. The relative standard deviations of B5-B6 are even larger than 40% in winter, as shown in Figure 10c,d. Besides the uncertainty of LSE, reflectance, and LST, it is also related to the radiation uncertainty caused by the coactions of the topographic uncertainty, and the low signal level in the spectral range of 1-5 µm. The TOA spectra of several typical areas and pixels, as shown in Figure 8, are compared with the onboard spectra gathered by the Gaofen-5 satellite. Figure 10 shows the TOA spectra comparisons among representative areas with different land covers. The ROI potions of cropland and water are marked in Figure 8b,d. The whole area of the rugged Mingsha mountain is taken as the ROI for the land cover of desert (sand). The mean spectra, with their adjacent interval within the one standard deviation range, are compared and shown in Figure 10a,c,d. One pixel of the road in the urban area is also selected as an example, representing the land cover of construction. In all four of the subfigures, the onboard and the simulated spectra have similar trends in spectral domain. The spectral differences between different imaging times are remarkable for all the four representative areas. In other words, the summer spectra and the winter spectra are distinguished within the one time of standard deviation areas. In addition, the cropland and the urban pixels have more vegetation covers in summer compared to winter. Accordingly, the simulated summer spectra in Figure 10a,b display obvious absorption features in the red spectral band, which is consistent with the onboard data. It is also obvious that the absolute value of the spectra standard deviation in the visual to SWIR bands are larger than the values in the MWIR and LWIR bands. The relative values of standard deviation in the spectral range of 1-5 μm are larger than other bands. The relative standard deviations of B5-B6 are even larger than 40% in winter, as shown in Figure 10c,d. Besides the uncertainty of LSE, reflectance, and LST, it is also related to the radiation uncertainty caused by the coactions of the topographic uncertainty, and the low signal level in the spectral range of 1-5 μm. The TOA spectral comparisons of three representative pixels in the Mingsha mountain areas are shown in Figure 11. The three pixels have the same reflectance spectra in the artificial scene, while suffering different degrees of surface relief. P1 is located near an urban area, with almost no topographic relief. The shade factor of P1 is equal to 1 for both of the two imaging times. P2 is located on the sunny side of a small slope. P3 has the largest slope angle compared with the other two points, and suffers shadow interference. The TOA spectral comparisons of three representative pixels in the Mingsha mountain areas are shown in Figure 11. The three pixels have the same reflectance spectra in the artificial scene, while suffering different degrees of surface relief. P1 is located near an urban area, with almost no topographic relief. The shade factor of P1 is equal to 1 for both of the two imaging times. P2 is located on the sunny side of a small slope. P3 has the largest slope angle compared with the other two points, and suffers shadow interference. In the winter scene, the shade factor of P3 is only 0.475, meaning that P3 receives just about half of the solar irradiance received by P1 and P2. The LST of P3 is 10 K lower than the LST of P1 and P2 in summer, while in the winter scene, the LST of P3 is about 30 K lower than the LST of P1 and P2. Consequently, the TOA spectra of P1 and P2 are obviously larger than those of P3. The TOA radiance simulation also realized reasonable modeling of the spectral difference between the two different imaging times. The TOA radiance spectra in the summer scene, with an imaging date of 14 July 2019, are larger than those in the winter scene on 21 December 2019. All the simulated spectra of the three points have a similar tendency in spectral domain to the corresponding onboard VIMS and AHSI data. Among all 12 bands of VIMS, the relative differences between the simulated and the onboard spectra of SWIR to MWIR spectral region are relatively larger than other spectral regions. This might be connected to calibration accuracy of the onboard data. In the winter scene, the shade factor of P3 is only 0.475, meaning that P3 receives just about half of the solar irradiance received by P1 and P2. The LST of P3 is 10 K lower than the LST of P1 and P2 in summer, while in the winter scene, the LST of P3 is about 30 K lower than the LST of P1 and P2. Consequently, the TOA spectra of P1 and P2 are obviously larger than those of P3. The TOA radiance simulation also realized reasonable modeling of the spectral difference between the two different imaging times. The TOA radiance spectra in the summer scene, with an imaging date of 14 July 2019, are larger than those in the winter scene on 21 December 2019. All the simulated spectra of the three points have a similar tendency in spectral domain to the corresponding onboard VIMS and AHSI data. Among all 12 bands of VIMS, the relative differences between the simulated and the onboard spectra of SWIR to MWIR spectral region are relatively larger than other spectral regions. This might be connected to calibration accuracy of the onboard data. To quantitatively compare the deviations between simulated and onboard data over rugged land surfaces, the RMSRDs of the Mingsha mountain pixels are calculated and shown in Table 4. It is obvious that the overall RMSRD of the summer scene is smaller than that of the winter scene. The smallest (largest) RMSRD is 3.35% (16.49%) of B12 (B8) in the summer scene, while the smallest (largest) RMSRD is 6.78% (37.31%) of B10 (B6) in the winter scene. The simulation accuracy of the proposed method is related to the  Table 4. It is obvious that the overall RMSRD of the summer scene is smaller than that of the winter scene. The smallest (largest) RMSRD is 3.35% (16.49%) of B12 (B8) in the summer scene, while the smallest (largest) RMSRD is 6.78% (37.31%) of B10 (B6) in the winter scene. The simulation accuracy of the proposed method is related to the imaging time, and varies with the wavelength. The RMSRDs of the near infrared to mid-wave infrared bands (NMWIR, B4 to B8) are relatively larger than other spectral regions for both of the two imaging times. A potential cause of the spectral and time variations of the deviations between simulated and onboard data is the radiometric calibration uncertainties of the VIMS. The restored TOA radiance between VIMS and AHSI are compared using the three typical ground points, shown in Figure 8. The original VIMS L1 data are restored to TOA radiance after the cross calibration. The restored TOA radiance spectra of AHSI are converted to equivalent spectra of the B1 to B6 bands using the SRFs of VIMS. The relative deviation between the VIMS and AHSI spectra are then calculated to represent the uncertainty of cross calibration over rugged areas, and are shown in Table 5. The results indicate that the relative deviations between sensors vary with imaging time and ground location. Therein, the relative deviations of the shaded pixel (P3) in B5 to B6 are larger than 150%. It is reasonable to believe that the winter scene of VIMS suffers larger uncertainties caused by cross calibration than the summer scene. Meanwhile, the bands of B4 to B6 suffer larger calibration uncertainties than other bands. Considering the statistical results in Tables 4 and 5, the simulation experiment indicates that the proposed simulation framework is reliable, because the RMSRD is less than 20% under most conditions, except for the NMWIR bands in winter with large solar zenith angles. In addition, the calibration uncertainty over a rugged surface is essential and warrants more attention in future quantitative applications.

Evaluation of Atmospheric and Topographic Contributions
To simulate the spectral imaging TOA data of rugged land surfaces under a given observing geometry condition in the visible to thermal domain, the proposed simulation framework introduces seven different radiative components, as defined in (2). Removing the influence of the upward direct transmittance, the atmospheric adjacency effect, and the path radiance, the BOA radiance is made up of five different components, as shown in (9). Both the TOA and BOA radiance are influenced by the atmosphere and terrain coupling effects. Although their mechanisms are different, the seven TOA components and the five BOA components are all variates of atmospheric and topographic parameters. Recently, it is a common cognizance that topographic parameters influence radiation signals, especially the downward solar and sky irradiance. However, there is still a lack of quantitative evaluations of components caused by different atmosphere and terrain coupling mechanisms over the visible to thermal domain, to the best of our knowledge. For the evaluation of atmospheric and topographic contributions on radiance signals, the various components in TOA and BOA radiance are calculated and compared. The relative contribution of components in the BOA and TOA radiance at the three typical pixels mentioned in the last section are expressed as percentage and shown in Figures 12 and 13 Recently, it is a common cognizance that topographic parameters influence radiation signals, especially the downward solar and sky irradiance. However, there is still a lack of quantitative evaluations of components caused by different atmosphere and terrain coupling mechanisms over the visible to thermal domain, to the best of our knowledge. For the evaluation of atmospheric and topographic contributions on radiance signals, the various components in TOA and BOA radiance are calculated and compared. The relative contribution of components in the BOA and TOA radiance at the three typical pixels mentioned in the last section are expressed as percentage and shown in Figure 12 and Figure 13, respectively. Thus, the component contribution order in the radiance signal of each pixel is clear.   Figure 12 illustrates that the reflection of the direct solar irradiance of the target is always the most significant component for the BOA radiance of the illuminated pixels (P1 and P2) in solar-reflected spectral region. The target emission is the most significant component in the spectral range from 4 μm to 12 μm, while the reflection of sky irradiance plays an important role through all the spectral bands, from visual to long-wave infrared (LWIR), especially for the shaded pixels, as shown in Figure 12f. The topographic irradiance effect is quite weak for the illuminated pixels, but remarkable for the shaded pixel. For the B6 of VIMS, with a central wavelength of 3.68 μm in the scene with the imaging date of 21 December 2019, the topographic irradiance contributes nearly 40% of the BOA radiance. Meanwhile, the trapping effect, which decreases when the wavelength increases, is more significant for the shaded pixel compared to illuminated pixels.
As shown in Figure 13, the atmospheric path radiance is the third significant component in TOA radiance, besides the target emission and the reflection of direct solar irradiance. The contribution of atmospheric path radiance is up to 76% of the condition at the first band of Figure 11f, which explains why the target recognition is harder in shadow areas. The adjacency effect on TOA radiance is the summation of the topographic irradiance, the trapping effect, and the atmospheric adjacency effect. Among the three adjacency effects, atmospheric adjacency effect is remarkable for all six of the conditions in Figure 13, while the topographic irradiance and the trapping effect contribute no more than 1% for the illuminated P1 and P2. As the contribution of topographic irradiance effect increases with the terrain shade, it constitutes 26% of the TOA radiance at B6 in Figure  13f. Although smaller contributions, caused by topographic irradiance and trapping effect, are observed in thermal domain than in the solar-reflected spectral domain, they are actually important for the shaded pixel P3, with a total contribution of 3% for B7.
The results in Figures 12 and 13 illustrate that the adjacency effects caused by atmosphere and terrain have small, but unstable, contributions to the BOA and TOA  Figure 12 illustrates that the reflection of the direct solar irradiance of the target is always the most significant component for the BOA radiance of the illuminated pixels (P1 and P2) in solar-reflected spectral region. The target emission is the most significant component in the spectral range from 4 µm to 12 µm, while the reflection of sky irradiance plays an important role through all the spectral bands, from visual to long-wave infrared (LWIR), especially for the shaded pixels, as shown in Figure 12f. The topographic irradiance effect is quite weak for the illuminated pixels, but remarkable for the shaded pixel. For the B6 of VIMS, with a central wavelength of 3.68 µm in the scene with the imaging date of 21 December 2019, the topographic irradiance contributes nearly 40% of the BOA radiance. Meanwhile, the trapping effect, which decreases when the wavelength increases, is more significant for the shaded pixel compared to illuminated pixels.
As shown in Figure 13, the atmospheric path radiance is the third significant component in TOA radiance, besides the target emission and the reflection of direct solar irradiance. The contribution of atmospheric path radiance is up to 76% of the condition at the first band of Figure 11f, which explains why the target recognition is harder in shadow areas. The adjacency effect on TOA radiance is the summation of the topographic irradiance, the trapping effect, and the atmospheric adjacency effect. Among the three adjacency effects, atmospheric adjacency effect is remarkable for all six of the conditions in Figure 13, while the topographic irradiance and the trapping effect contribute no more than 1% for the illuminated P1 and P2. As the contribution of topographic irradiance effect increases with the terrain shade, it constitutes 26% of the TOA radiance at B6 in Figure 13f. Although smaller contributions, caused by topographic irradiance and trapping effect, are observed in thermal domain than in the solar-reflected spectral domain, they are actually important for the shaded pixel P3, with a total contribution of 3% for B7.
The results in Figures 12 and 13 illustrate that the adjacency effects caused by atmosphere and terrain have small, but unstable, contributions to the BOA and TOA radiance signals. Whether these adjacency effects could be ignored when using a specified imaging system becomes an important problem in the quantitative remote sensing. Taking the noise level assessment of VIMS as an example, the relative contribution of adjacency effects on TOA signal of different pixels are compared with the multiplicative inverse of SNR, as shown in Figure 14. The NE∆T values for B6 to B12 in Table 3 are converted to SNR at a signal level, defined by blackbody radiance with a temperature of 300 K. The multiplicative inverse of SNR is regarded as the relative radiation resolution of the imager. Figure 14 shows that the total relative contributions of the three adjacency effects on the TOA radiance are larger than the multiplicative inverse of SNR, for most spectral regions. The gaps between the relative contributions and the multiplicative inverse of SNR are obviously larger in the solar-reflected spectral domain than those in the thermal domain. The logarithmic axis method is used for the y-axis of this figure to display a more detailed contrast between plots in the thermal spectral domain. The spectral relative contributions of the three adjacency effects for the three typical pixels all have their smallest values at the band of B9, with a central wavelength of 8.19 µm. Only the contributions of the illuminated P1 and P2 pixels at B9 are smaller, or approximate to, the multiplicative inverse of SNR. In other words, only for the illuminated P1 and P2 pixels at B9 can the adjacency effects be ignored in the measured TOA signals. The relative contributions of the three adjacency effects are smaller than 5% over the thermal domain. It implies that these effects are not accurately measured when the restored TOA data has a relative measurement error of 5%, leading to uncertainty in remote sensing products.
Remote Sens. 2022, 14, x FOR PEER REVIEW 24 of 30 radiance signals. Whether these adjacency effects could be ignored when using a specified imaging system becomes an important problem in the quantitative remote sensing. Taking the noise level assessment of VIMS as an example, the relative contribution of adjacency effects on TOA signal of different pixels are compared with the multiplicative inverse of SNR, as shown in Figure 14. The NE△T values for B6 to B12 in Table 3 are converted to SNR at a signal level, defined by blackbody radiance with a temperature of 300 K. The multiplicative inverse of SNR is regarded as the relative radiation resolution of the imager. Figure 14 shows that the total relative contributions of the three adjacency effects on the TOA radiance are larger than the multiplicative inverse of SNR, for most spectral regions. The gaps between the relative contributions and the multiplicative inverse of SNR are obviously larger in the solar-reflected spectral domain than those in the thermal domain. The logarithmic axis method is used for the y-axis of this figure to display a more detailed contrast between plots in the thermal spectral domain. The spectral relative contributions of the three adjacency effects for the three typical pixels all have their smallest values at the band of B9, with a central wavelength of 8.19 μm. Only the contributions of the illuminated P1 and P2 pixels at B9 are smaller, or approximate to, the multiplicative inverse of SNR. In other words, only for the illuminated P1 and P2 pixels at B9 can the adjacency effects be ignored in the measured TOA signals. The relative contributions of the three adjacency effects are smaller than 5% over the thermal domain. It implies that these effects are not accurately measured when the restored TOA data has a relative measurement error of 5%, leading to uncertainty in remote sensing products.

Figure 14.
Overall contributions of atmospheric and topographic adjacency effects on TOA radiance.
As a preview of the errors on TOA brightness temperature products caused by ignoring the adjacency effects, the deviation of TOA brightness temperature (DTBT) of simulated images with, and without, the related three components are shown in Figure  15. Obviously, the influence of adjacency effects on DTBT is spectrally and spatially variable. The images show that the DTBT of the pixels in the valleys, with small sky view factors, is relatively higher than other pixels. In the flat and peak area, the DBT is mainly distributed in a range of 0.2 to 0.7 K for all four LWIR bands. For the valley pixels, the DTBT is larger than 1 K. The largest DTBT of 1.75 K appears at B10 on the day of 14 July 2019. Although the simulations are just taken under clear and sunny conditions, they also have referential values for people to recognize the uncertainty associated with adjacency effects over rugged areas in the visible to thermal spectral domain. As a preview of the errors on TOA brightness temperature products caused by ignoring the adjacency effects, the deviation of TOA brightness temperature (DTBT) of simulated images with, and without, the related three components are shown in Figure 15. Obviously, the influence of adjacency effects on DTBT is spectrally and spatially variable. The images show that the DTBT of the pixels in the valleys, with small sky view factors, is relatively higher than other pixels. In the flat and peak area, the DBT is mainly distributed in a range of 0.2 to 0.7 K for all four LWIR bands. For the valley pixels, the DTBT is larger than 1 K. The largest DTBT of 1.75 K appears at B10 on the day of 14 July 2019. Although the simulations are just taken under clear and sunny conditions, they also have referential values for people to recognize the uncertainty associated with adjacency effects over rugged areas in the visible to thermal spectral domain.

Discussion
Although the validation results demonstrate that the proposed atmosphere and terrain coupling simulation framework over visible to thermal domain could gather consistent and accurate simulations of onboard images, there are some issues that require attention when dealing with simulation and data processing tasks.

Topographic Error Accumulation
In the proposed simulation framework, all the TOA and BOA radiance components, and the simulated LST, are related to topographic features. The slope angle, aspect angle, cosine of SIA, shade factor, sky view factor, terrain view factor, the cosine of the viewing angle (VA, angle between the normal of pixel and the viewing direction of sensor), and the LST are collectively called as topographic variable factors here. Errors of all these topographic variable factors would accumulate in the simulated BOA and TOA spectral radiance. The terrain view factor in (25), with a value from 0 to 1, is introduced to depict the degree of adjacency illumination for discussion in this section.

cos 2 t sky
In order to identify the major topographic error sources, the correlation coefficients between the input topographic variable factors and the absolute values of relative deviation of the simulated and onboard images are shown in Figure 16. As the input atmospheric coefficients are spatially interpolated, according to the ground elevation, the DEM data is also discussed as a variable factor here. On the whole, the value of the topographic correlation of the summer scene is greater than that of the winter scene, which may be related to the different calibration accuracies of onboard data. However, Figure 16 provides a coherent view of the influence rank on overall accuracy of the proposed model over the visible to thermal domain. The influence rank, in descending order, is cosine of SIA, shade factor, slope, LST, cosine of VA, sky view factor, and terrain view factor, while the influence on simulation accuracy of the DEM and aspect angle are relatively nonlinear and indirectly transferred by other topographic variable factors. These key topographic

Discussion
Although the validation results demonstrate that the proposed atmosphere and terrain coupling simulation framework over visible to thermal domain could gather consistent and accurate simulations of onboard images, there are some issues that require attention when dealing with simulation and data processing tasks.

Topographic Error Accumulation
In the proposed simulation framework, all the TOA and BOA radiance components, and the simulated LST, are related to topographic features. The slope angle, aspect angle, cosine of SIA, shade factor, sky view factor, terrain view factor, the cosine of the viewing angle (VA, angle between the normal of pixel and the viewing direction of sensor), and the LST are collectively called as topographic variable factors here. Errors of all these topographic variable factors would accumulate in the simulated BOA and TOA spectral radiance. The terrain view factor in (25), with a value from 0 to 1, is introduced to depict the degree of adjacency illumination for discussion in this section.
In order to identify the major topographic error sources, the correlation coefficients between the input topographic variable factors and the absolute values of relative deviation of the simulated and onboard images are shown in Figure 16. As the input atmospheric coefficients are spatially interpolated, according to the ground elevation, the DEM data is also discussed as a variable factor here. On the whole, the value of the topographic correlation of the summer scene is greater than that of the winter scene, which may be related to the different calibration accuracies of onboard data. However, Figure 16 provides a coherent view of the influence rank on overall accuracy of the proposed model over the visible to thermal domain. The influence rank, in descending order, is cosine of SIA, shade factor, slope, LST, cosine of VA, sky view factor, and terrain view factor, while the influence on simulation accuracy of the DEM and aspect angle are relatively nonlinear and indirectly transferred by other topographic variable factors. These key topographic variable factors are extracted from the DEM data. Slope and aspect angles are calculated directly from DEM, and they are mathematically correlative to each other. There would be error accumulation of topographic variable factors when radiative factors over rugged areas are calculated with these factors. Therefore, the quantitative description of errors and correlations of topographic variable factors is an interesting and important task in future studies to explain the BOA product uncertainties. variable factors are extracted from the DEM data. Slope and aspect angles are calculated directly from DEM, and they are mathematically correlative to each other. There would be error accumulation of topographic variable factors when radiative factors over rugged areas are calculated with these factors. Therefore, the quantitative description of errors and correlations of topographic variable factors is an interesting and important task in future studies to explain the BOA product uncertainties.

The Angle and Scale Effects
This study is built on the DEM data that describes the surface terrain and the 1D plane-parallel atmosphere assumption, which segments the atmospheric medium into horizontal layers. The validation experiments were taken from the simulation of the imaging data of a satellite observing the Earth along nadir direction. The GSD of DEM data used in the experimental study in this paper is 12.5 m. The whole atmospheric medium of 100 km in altitude was segmented into 28 layers, automatically following the segmentation rule of the auxiliary profile data. Currently, the validation experiments have not been conducted using airborne data, with GSD less than 10 m. In this circumstance, the angle and scale effects become a practical problem. A higher resolution of DEM data is required to solve the problems caused by composite slope effects and to gather accurate topographic parameters.
It is worth noting that the atmospheric PSF, as shown in (16), is a function of the viewing angle. For the nadir observation over flat surfaces, the atmospheric PSF is symmetrical with the highest value near the target pixel, because the aerosol scattering phase function peaks strongly in the forward direction. It is believed that the shape of atmospheric PSF is not important to the radiative results, as shown in previous literature such as [32]. Figure 17 presents the modeled atmospheric adjacency effect convolution function for typical pixel P3, under the nadir spaceborne viewing condition, and a tilt airborne viewing condition with view zenith angle of 20°. For the airborne condition, the integral calculation of atmospheric PSF was conducted with a refined vertical interval to 10 m, to decrease the extra error in numerical calculation. The shapes of these convolution functions are affected by local terrain and are not perfectly symmetrical. There is still a sharp peak near the target pixel, due to the nadir observation of spaceborne VIMS, as shown in Figure 17a. Nevertheless, for the airborne condition with tilt viewing directions and short transmission paths from ground to sensor, the geometric gravity center of the kernel function F moves away from the target location, as shown in Figure 17b. A singlepeaked hypothesis is invalid for atmospheric PSF for the tilt airborne viewing conditions. This phenomenon is related to angle and scale effects, and would be serious for the edge pixels in the sensor's field of view. Under the tilt viewing conditions, the ground-leaving

The Angle and Scale Effects
This study is built on the DEM data that describes the surface terrain and the 1D planeparallel atmosphere assumption, which segments the atmospheric medium into horizontal layers. The validation experiments were taken from the simulation of the imaging data of a satellite observing the Earth along nadir direction. The GSD of DEM data used in the experimental study in this paper is 12.5 m. The whole atmospheric medium of 100 km in altitude was segmented into 28 layers, automatically following the segmentation rule of the auxiliary profile data. Currently, the validation experiments have not been conducted using airborne data, with GSD less than 10 m. In this circumstance, the angle and scale effects become a practical problem. A higher resolution of DEM data is required to solve the problems caused by composite slope effects and to gather accurate topographic parameters.
It is worth noting that the atmospheric PSF, as shown in (16), is a function of the viewing angle. For the nadir observation over flat surfaces, the atmospheric PSF is symmetrical with the highest value near the target pixel, because the aerosol scattering phase function peaks strongly in the forward direction. It is believed that the shape of atmospheric PSF is not important to the radiative results, as shown in previous literature such as [32]. Figure 17 presents the modeled atmospheric adjacency effect convolution function for typical pixel P3, under the nadir spaceborne viewing condition, and a tilt airborne viewing condition with view zenith angle of 20 • . For the airborne condition, the integral calculation of atmospheric PSF was conducted with a refined vertical interval to 10 m, to decrease the extra error in numerical calculation. The shapes of these convolution functions are affected by local terrain and are not perfectly symmetrical. There is still a sharp peak near the target pixel, due to the nadir observation of spaceborne VIMS, as shown in Figure 17a. Nevertheless, for the airborne condition with tilt viewing directions and short transmission paths from ground to sensor, the geometric gravity center of the kernel function F moves away from the target location, as shown in Figure 17b. A single-peaked hypothesis is invalid for atmospheric PSF for the tilt airborne viewing conditions. This phenomenon is related to angle and scale effects, and would be serious for the edge pixels in the sensor's field of view. Under the tilt viewing conditions, the ground-leaving radiance of more adjacency pixels would be scattered into the sensor's IFOV. It is necessary to distinguish the local shape of atmospheric PSF varying with pixel locations. radiance of more adjacency pixels would be scattered into the sensor's IFOV. It is necessary to distinguish the local shape of atmospheric PSF varying with pixel locations. However, these potential angle and scale effects on airborne conditions could be modeled using the proposed simulation framework, with appropriate input data. The proposed simulation framework would benefit future improvement of atmospheric correction methods over rugged land surfaces.

Challenges for Retrieval Applications
The forward simulation method is proven effective, as the results agree reasonably well with satellite observations statistically over a give rugged land scene. However, the proposed model could not be directly used for backward retrieval of LSE, reflectance, and LST. There are still some challenges needing attention before developing the atmospherictopographic coupling correction methods.
(1) The measurement function of TOA radiance is an ill-posed function, with more independent variables than dependent variables. The backward retrieval needs to be simplified to the normal functions, or solved using iteration. Thus, parameterization of the initial LSE, reflectance, and LST with reliable prior knowledge, is the first challenge.
(2) The proposed model describes the 3D geometry of land surface using the public DEM data, with a rough GSD of 10 m to 100 m. Accurate descriptions of the terrain futures of pixels in the retrieval tasks requires the GSD of imaging data larger than the GSD of DEM data. Thus, the generation of fine and accurate geometry parameters is another challenge when dealing with complex 3D surfaces, such as the canopy and the buildings.
(3) The atmospheric-topographic coupling correction in the scene with clouds is also a large challenge. This study is conducted without the 3D effect of clouds, which are essential radiation sources for the land surface, as well as obstacles for both upward and downward radiation transfer. The 1D cloudy atmosphere influence on the atmospheric flux and the transmission could be modelled using the MODTRAN module in the proposed simulation framework. Moreover, the priori label layers of water clouds, cirrus clouds, and their cast shadow over the rugged surface are necessary in the retrieval tasks.

Conclusions
As the continuation of our previous work of the earth-atmosphere coupling model for solar-reflected hyperspectral systems [29], the model was improved and expanded to the entire spectral range, from visible to thermal. The atmosphere and terrain coupling simulation framework for high-resolution spectral imaging of heterogeneous land However, these potential angle and scale effects on airborne conditions could be modeled using the proposed simulation framework, with appropriate input data. The proposed simulation framework would benefit future improvement of atmospheric correction methods over rugged land surfaces.

Challenges for Retrieval Applications
The forward simulation method is proven effective, as the results agree reasonably well with satellite observations statistically over a give rugged land scene. However, the proposed model could not be directly used for backward retrieval of LSE, reflectance, and LST. There are still some challenges needing attention before developing the atmospherictopographic coupling correction methods.
(1) The measurement function of TOA radiance is an ill-posed function, with more independent variables than dependent variables. The backward retrieval needs to be simplified to the normal functions, or solved using iteration. Thus, parameterization of the initial LSE, reflectance, and LST with reliable prior knowledge, is the first challenge.
(2) The proposed model describes the 3D geometry of land surface using the public DEM data, with a rough GSD of 10 m to 100 m. Accurate descriptions of the terrain futures of pixels in the retrieval tasks requires the GSD of imaging data larger than the GSD of DEM data. Thus, the generation of fine and accurate geometry parameters is another challenge when dealing with complex 3D surfaces, such as the canopy and the buildings.
(3) The atmospheric-topographic coupling correction in the scene with clouds is also a large challenge. This study is conducted without the 3D effect of clouds, which are essential radiation sources for the land surface, as well as obstacles for both upward and downward radiation transfer. The 1D cloudy atmosphere influence on the atmospheric flux and the transmission could be modelled using the MODTRAN module in the proposed simulation framework. Moreover, the priori label layers of water clouds, cirrus clouds, and their cast shadow over the rugged surface are necessary in the retrieval tasks.

Conclusions
As the continuation of our previous work of the earth-atmosphere coupling model for solar-reflected hyperspectral systems [29], the model was improved and expanded to the entire spectral range, from visible to thermal. The atmosphere and terrain coupling simulation framework for high-resolution spectral imaging of heterogeneous land surfaces over visible to thermal domain was established in this study. The framework, in combination with the heat balance model of land surfaces, realized the simulation of high-resolution imaging signals including the LST, the BOA radiance, and the TOA radiance over rugged areas. There are two advantages of the simulation methodology: (I) it analytically models the topographic irradiance, trapping effect, and atmospheric adjacency effect with three spatialspectral convolution kernels related to both atmospheric and terrain features, avoiding the commonly used assumption of independence between atmospheric and topographic adjacency effects; (II) it provides a novel method to establish layer atmosphere coefficients look-up tables based on multiple output files of MODTRAN, which benefits the imaging simulation for various viewing altitude and angles.
The simulation framework, and the produced data, could be used for the performance predictions of earth-observing systems, and serve to optimize spectral imaging missions. The accuracy of the TOA radiance simulation was validated using China's Gaofen-5 VIMS data. The simulation experiment indicates that the proposed simulation framework is reliable. The spectral root mean square of relative deviations between simulated radiance and onboard radiance is less (larger) than 17% (3%) and 38% (6%) in the summer and winter scene, respectively. Some statistical evaluations of the atmospheric and topographic effects were also implemented. It indicates that the uncertainty associated with atmosphere and terrain coupling effects should be quantified under the whole cycle, from the concept design stage to the application stage, for a spectral imaging mission of earth.