MRRT: Modeling of Multiple Reflections of Radiation Between Terrains Applied on the Lunar Surface Region

The lunar surface has a stable luminosity. To use the Moon as a calibration standard, the Robotic Lunar Observatory (ROLO) program models the integration of the radiance of the entire lunar surface. However, the albedos of the mare and the highlands are very different. The modeling based on the lunar global irradiance/reflected radiance is bound to result in higher uncertainty. In contrast, if the local calibration of the lunar surface is adopted, the lunar complex topography effect cannot be ignored. This article presents a new model for quantifying multiple reflections of radiation between terrains (MRRT). The relationship between the bidirectional reflectance factor (BRF) of the observed pixel and the true microtopography reflectance is established, which shows that the BRF is mainly influenced by the true topography reflectance, the terrain undulation, the incident irradiance on the topography surface, and the masking in the observation direction. The new model applied on the lunar surface obtains clearer terrain details. The inversion reflectance of the Chang’e-3 landing area is closer to the reflectance measured in situ, and the reflectance curves of the Apollo 16 landing area are almost consistent under different illumination observation geometries. This shows that the MRRT model can effectively eliminate the topographic effect. Compared with the ROLO model, the MRRT model does not restrict the specific selection, so it can select a region with a uniform material distribution, small albedo difference, and low topography undulation to establish the lunar surface radiometric calibration field with the advantage of providing stable radiation characteristics.

MRRT: Modeling of multiple reflections of radiation between terrains applied on the lunar surface region Yunfei Liu, Qiang Guo and Guifu Wang Abstract-The lunar surface has a stable luminosity.To use the Moon as a calibration standard, the Robotic Lunar Observatory (ROLO) program models the integration of the radiance of the entire lunar surface.However, the albedos of the mare and the highlands are very different.The modeling based on the lunar global irradiance/reflected radiance is bound to result in higher uncertainty.In contrast, if the local calibration of the lunar surface is adopted, the lunar complex topography effect cannot be ignored.This paper presents a new model for quantifying multiple reflections of radiation between terrains (MRRT).The relationship between the bidirectional reflectance factor (BRF) of the observed pixel and the true microtopography reflectance is established, which shows that the BRF is mainly influenced by the true topography reflectance, the terrain undulation, the incident irradiance on the topography surface, and the masking in the observation direction.The new model applied on the lunar surface obtains clearer terrain details.The inversion reflectance of the Chang'e-3 landing area is closer to the reflectance measured in situ, and the reflectance curves of the Apollo 16 landing area are almost consistent under different illumination observation geometries.This shows that the MRRT model can effectively eliminate the topographic effect.Compared with the ROLO model, the MRRT model does not restrict the specific selection, so it can select a region with a uniform material distribution, small albedo difference, and low topography undulation to establish the lunar surface radiometric calibration field with the advantage of providing stable radiation characteristics.
Index Terms-Effect of microtopography inside remote sensing observation pixels, multiple-reflection-process modeling of radiation between terrains, and inversion of lunar surface reflectance

I. INTRODUCTION
T HE Moon is a nearby atmosphere-less body that reflects solar radiation [1].For fixed illumination and fixed observation geometry, the Moon can be considered photometrically stable up to 10 −8 per annum in irradiance [2].Therefore, the Moon can serve as an ideal external calibration source for remote sensing instruments in the reflected solar band.The United States Geological Survey ran the Robotic Lunar Observatory (ROLO) ground-based observational program to establish an empirical lunar radiation model, which is based on observational data and modeling the integration of the radiance of the entire lunar surface at different times [3].Although the ROLO model is currently the most accurate lunar radiation model worldwide, it still has an uncertainty of 5% to 10% This work was supported by National Natural Science Foundation of China under Grant 41875037 and project for FY-4 ground-based visible and infrared lunar observation system of Shanghai Institute of Technical Physics, Chinese Academy of Sciences (No. O9KCE033N3).(Corresponding author: Qiang Guo.) Manuscript received April 19, 2021; revised August 16, 2021 [4].One of the reasons is that the ROLO model reflectance is adjusted by using the laboratory reflectance of returned Apollo samples.The reflectance of returned Apollo soil samples in the laboratory is significantly different from that from lunar observations because of the difference in gravity between the Moon and Earth [5].To improve the ROLO model, Zhang et al. [5] proposed using the mean equigonal albedo to replace the reflectance of Apollo soil samples; the new model's irradiance values in the visible and near-infrared (VNIR) bands are closer than those of the traditional ROLO model to the observations from the Moderate Resolution Imaging Spectroradiometer (MODIS) and Sea-viewing Wide Field-ofview Sensor (SeaWiFS).Sun et al. [6] developed a new model based on MODIS instrument observations to compensate for the shortcomings of the ROLO model in lunar irradiance measurement, which significantly improved the MODIS lunar calibration results of the entire mission.In fact, the Moon is not an ideal Lambertian body, and the albedos of different regions of the lunar surface are quite different (see Table I [7]).
A modeling method based on the global irradiance/reflected radiance of the Moon is bound to result in higher uncertainty.Large and small impact craters are densely distributed on the lunar surface [8].At the meter to hectometer scales, there are obvious differences in the median bidirectional slope, root-mean-square (RMS) height, and median absolute slope between mares and highlands [9].The details of the luminosity properties of the lunar surface are complicated by the macroscopic roughness of the Moon [10].Rugged terrain often alters illumination and viewing geometry and generates a relief shadow, observation masking, and multiple scattering, which result in intense topographic dependence on the total incident reflectance or radiance [11]- [13].According to different terrane types of the lunar surface, Wu et al. [14] divided the moon into four classes of albedo types and established a lunar irradiance model based on Chang'e-1 imaging interferometer (IIM) data; however, their model has some shortcomings in that the wavelength range is too narrow, and the model errors in high-latitude areas and border areas are relatively large.If a local area on the lunar surface is used for calibration, it is necessary to consider incorporating topographic data to study the problem of rim modeling and to improve the irradiance simulation accuracy [14].
At present, empirical models such as the Lommel-Seeliger model [15] and Sandmeier model [16] are often adopted for the calibration of lunar photometric observations.Alternatively, the Hapke radiative transfer model with simplified parameters or other models can be used [17].Empirical models are often prone to overcorrection due to their simple parameters [18], [19], but the Hakpe model [20]- [23] is difficult to apply due to its numerous parameters, mathematical coupling, and high requirements regarding observation data [17].With the improvement of computing power, 3D computer simulation models play an increasingly important role in the study of the radiation characteristics of complex surfaces [24].The main principles of computer modeling include ray tracing and radiosity methods.Using statistical analysis by computers, the distribution of lunar soil particles can be simulated.Muinonen et al. [25] used the Monte Carlo ray tracing method to analyze the single-scattering albedos and phase functions, local surface roughness, and regolith porosity of specific lunar mare regions imaged by the Advanced Moon micro-Imager Experiment (AMIE) camera onboard the ESA SMART-1 mission.Wong et al. [26] combined Monte Carlo ray tracing and Hapke models to model reflectance considering both (surface) large scale and (interparticle) microscale effects.However, Wong used a typical simple illumination model, the Phong model, which does not consider reflection between terrain surfaces.
In ground remote sensing detection, Pory et al. [27] assumed that the surface is Lambertian and proposed that the radiation received by a surface consists of direct solar radiation, scattered atmospheric radiation, and adjacent terrain-reflected radiation.The reflected radiation from the surrounding terrain is defined as the sum of the solar radiation reflected on the target pixel by other visible pixels (see Equation ( 4)).Nevertheless, the calculation of reflected radiation from the surrounding terrain only considers the first-order scattering effect between terrains.The attenuation of the signal between two adjacent slopes and second-order reflection are neglected [27].In remote sensing on Earth, the multiple reflections of radiation between terrains are often attenuated greatly due to the influence of air, but for the lunar surface with atmosphereless, the multiple reflections of radiation between terrains should be given more attention.
In this paper, a new model for quantifying multiple reflections of radiation between terrains (MRRT) is presented.Based on the adjacent terrain irradiance formula of the first-order reflection proposed by Proy, the second-order to the nth-order reflections of radiation between terrains are derived.Moreover, the relationship between the bidirectional reflectance factor (BRF) of the observed pixel and the true microtopography reflectance is established, which shows that the BRF is mainly influenced by the true topography reflectance, the terrain undulation, the incident irradiance on the topography surface, and the masking in the observation direction.
In Section II, the modeling process of multiple reflections of radiation between terrains is described in detail.In Section III, the experimental area and data selected in this paper are introduced.In Section IV, we apply the new model to the experimental area and compare it with the dataset to discuss its results.Finally, in Section V, we summarize the content of this paper and prospect the application scenarios of the new model.

II. THEORY
Direct solar radiation over rugged terrain is the most important component of the total incident radiation that illustrates the surface [13].In addition, the adjacent terrain-reflected irradiance increases the total radiation reaching the slope surface [27].As shown in Fig. 1, M and P are mutually visible.When point P has a reflected radiance greater than zero, point M must be reflected by point P ; that is, M receives both incident radiation from the sun and radiation reflected from point P .Obviously, for the lunar surface with almost no atmospheric attenuation, the reflection of radiation between terrains (such as points P and M ) is performed multiple times.
Assuming that the total solar incident radiation flux of the target terrain is Φ sun (limited area size and no other incident radiation source), the total radiation flux of the first reflection is Φ 1 , and after k reflections, the total radiation flux is Φ k .The relationship among Φ sun , Φ 1 and Φ k can be expressed by the following formula: and where Φ k infinitely approaches zero.The target terrain is divided into equal space intervals so that multiple microareas with different slope, aspect, and elevation values are formed.This has the advantage of simplifying the shape of the terrain to discuss how the radiation varies from  microarea to microarea.As shown in Fig. 2, every reflection of radiation between terrains will cause changes in incident radiation on the surface of microareas.We call this change in incident radiation on the surface of the microarea an "update" of incident radiation on the microarea surface.
A. Updating the incident radiation on the surface of a microarea 1) One-time reflection between terrains: The surfaces with different reflectance are assumed to be Lambertian.As shown in Fig. 1, the radiance received by point M from point P can be written as [27]: where dS M and dS P are the areas of pixels M and P , respectively; T M and T P are the angles between the normal to the ground and the line M P ; L P is the luminance of P ; and r M P is the distance between M and P .If ρ P is the reflectance of pixel P illuminated by the irradiance E P , then L P = ρ P E P /π.The subscript P → M represents the transfer of radiation from P to slope M .
Therefore, the total irradiance received by slope M from all "visible" P [27] can be written as: 2) Multiple reflections between terrains: As noted in the previous section, if ρ M Φ M is not 0 (where ρ M is the reflectance of M and Φ M is the total incident radiation flux received by the slope M from other visible slopes after onetime reflection), then ρΦ M will continue to participate in the next reflection with other slopes.Therefore, the incident irradiance on each microarea surface after each radiation reflection can be obtained as follows: where E M (n) represents the total incident irradiance received by the slope M surface for the nth time. Let where Γ M Pj can be considered the visible radiation factor between P j and M .Obviously, Γ M Pj does not vary with the number of reflections.Equation (5b) can be simplified as: Then, after the multireflection between terrains, the total reflected radiance L M ref of slope M to the sky can be expressed as: Equation (8) shows that the solar direct incident irradiance on the slope surface, the mutual visibility between slopes and the reflectivity of slopes directly affect the multiple reflection process and results.

B. Direct solar incident irradiance on microtopography surfaces
Changes in altitude affect the distribution of solar radiation, resulting in sunlit and shaded areas that correspond to bright and dark pixels in remote sensing images [28].In general, on slopes of rugged terrain, the most important variable controlling incident radiation is the local solar illumination angle [29].If the sun is not hidden by a local horizon, the local illumination angle θ s on a slope S with azimuth A is given by: where θ 0 is the illumination angle on a horizontal surface and ϕ 0 is the azimuth of illumination.
Shadows refer to regions lacking direct solar illumination, which can be attributed to two reasons: opposing the sun and cast by obstructions [30].As shown in Fig. 3, for any slope, the received solar irradiation can be divided into three types: (i) fully irradiated, (ii) partially irradiated, and (iii) no irradiation.First, it is determined whether the two ends of a slope can be irradiated, and the result is used to establish the actual irradiated condition of the slope.
1) Incident irradiance of a one-dimensional terrain surface: We form grid points with elevation values.For a grid row, an elevation function z is defined for the points j = 0, 1, . . ., N − 1.Since the points are evenly spaced, the abscissa is specified by j∆h [29].A binary factor Θ is established to indicate whether elevation point z j is irradiated.When Θ = 1, z j is irradiated; when Θ = 0, z j is not irradiated (see Fig. 3).
Therefore, for all 0 ≤ i < N , when the direction of solar incidence is the opposite of the increasing direction of i: for all i < j < N , if where θ h is the solar elevation angle.There is at least one point j < N that can be connected to i to form a new slope so that the elevation angle of the new slope is greater than that of the sun.This indicates that point i cannot be irradiated; thus, For all 0 ≤ i < i + 1 < N , let dS i denote the slope consisting of the ith and the (i+1)th elevation values.If Θ i = 0, dS i is not irradiated at all; if Θ i = 1 and Θ i+1 = 1, dS i is fully irradiated.Then, the irradiance formula can be expressed as: If Θ i = 1, Θ i+1 = 0, then dS i is partially irradiated, as shown in Fig. 3.The size of the irradiated area depends on the maximum occlusion point corresponding to z i+1 (denoted as z k ).The critical point at which the slope dS i receives solar radiation can be obtained from the intersection of the line and z i z i+1 of solar rays passing through the point z k .Assuming that the critical point is z d , the irradiance formula for this irradiated area is: where E sun is the direct solar irradiance, ⃗ N is the normal to the terrain and ⃗ S is the solar angle.

2) Incident irradiance of a two-dimensional terrain surface:
The topographic effects on solar irradiance are mainly variations in illumination angle and shadowing from local horizons [29].Due to the irregular nature of the ground in rugged areas, the sky dome overlying a surface is not the integrated hemisphere of a horizontal surface [31].Dozier [32] proposed determining the local horizon information from a grid.At any location, the portion of the overlying hemisphere that is obscured by terrain is: where h[θ] is the horizon angle in the direction θ.These horizons, however, are difficult to compute because, unlike slope and azimuth, they cannot be generated from information restricted to the immediate neighborhood of a point.
By rotating a grid in direction ϕ 0 , we reduce the horizon problem to its one-dimensional equivalent.Our interest is the angle to the horizon from any point in any direction, but we formulate the problem by determining the coordinates of the points that form the horizons.Minor errors in the elevation grid can therefore shift the "answer", i.e., the coordinates of the horizon point, by a considerable distance, but minor errors do not cause much variation in the end result, which is the angle to the horizon.When available digital terrain grids are sufficiently smooth, interpolation does not change the results [32].In Appendix A, we discuss the effect of different interpolation methods on the error caused by image rotation.
A suitable interpolation method is selected for the rotation of the grid so that the relative azimuth of the incident solar  ray is 0 • .In this paper, the bicubic interpolation method with the minimum error after the two times rotation is used.The kernel function of bicubic interpolation [33] is as follows: where a is equal to −0.5.The rotated grid is calculated using the one-dimensional terrain surface incident irradiance method, and then the calculated matrix is rotated in the inverse direction of ϕ 0 to obtain the irradiance of the original grid.
It is assumed that the final solar direct incident irradiance is represented by matrix E 0 , which can be written as: E rc represents the actual solar incident irradiance received by the grid points located in row r and column c.

C. Visible radiation factor between terrains
To quantify the magnitude of the new incident radiation on the surface of a microarea after each reflection, we need to determine in advance what other surfaces can reflect the target.In this paper, the mutual reflection between surfaces is reduced to the mutual visibility between grid points.One grid point represents a flat cylinder whose height is the elevation value of the terrain, as shown in Fig. 4.
For a given terrain, the mutual reflection relationship between any slope pair is also certain.We first calculate whether there is a reflection between any two slopes; the calculation result can be stored in advance so that it can be called directly in the later operation.The binary matrix F (r,c) stores the reflection relations between the grid point at position (r, c) and all other slope points, and the matrix elements are 0 or 1.As shown in Fig. 5, F (1,1) is used as an example: when F (1,1) (1, 1) = 0, the DEM pixel at position (1, 1) in the grid matrix includes no self-reflection; when F (1,1) (1, 2) = 1, reflection occurs between DEM pixels at position (1, 2) and pixels at position (1, 1).
We fit the elevation value of the "position" on the linear path between the "radiation source point" and the "target receiving point" by the adjacent point values and determine whether the radiation rays are "blocked".As shown in Fig. 4(b), there are not always true elevation points on the radiating paths of A → B and A → C (i.e., there are no elevation points at the double circle).We define the position of the elevation value passing on the radiation path between the radiation source and the target receiving point as the "intermediate slope point", which may or may not be present in the DEM.The number of intermediate slope points depends on the maximum distance between the radiation source and the target receiver in the row and col directions.
As shown in Fig. 6, the Oxyz spatial Cartesian coordinate system is established, where x and y correspond to the row is the elevation value of (x d , y2).Then, the formula for z d is: For y a < y d < y b , x d is the intersection of the perpendicular of (0, y d ) in the direction of x and the projection of AB in the two-dimensional plane of Oxy.x1 = f loot(x d ), and x2 = ceil(x d ), where z1 is the elevation value of (x 1 , y d ) and z2 is the elevation value of (x2, y d ).Then, the formula for z d is: Therefore, it is possible to calculate whether any point D on the "radiation path" of AB occludes AB.The judgment process is similar to that expressed in Equation (10).Alternatively, this judgement can be based on the position of z d (i.e., whether z d is above or below line AB).
Based on the binary matrix F of each pixel position, the visible radiation factor Γ of terrain corresponding to each pixel position can be obtained by combining Equation ( 6).

D. Matrix representation of multiple reflections between terrains
According to the content in the previous three sections, the multiple-reflection process between terrains is expressed in the form of a matrix.The microarea reflectance is assumed to be ρ, and Lambert reflection is considered.Let the matrix T n store the incident irradiance on the surface of each DEM pixel of the nth-order, as shown in Equations ( 19) and ( 20): T T n−1 ⊙ Γ (r,c) represents the multiplication of each element in the corresponding position of matrix T n−1 and matrix Γ(r, c), resulting in a matrix of the same size.sum(T n−1 ⊙ Γ (r,c) ) denotes adding all the elements of the resulting matrix.Each element in T n gradually tends to 0 over multiple radiative reflections.Then, the final total reflected radiance of the whole DEM image to the sky L ref can be expressed as: E. The total reflected radiance in the observed direction Assume that the illumination of each microarea at the incident angle is given by matrix B s and that the visibility of each microarea at the observation angle is given by matrix B v .According to the principle of reciprocity of angles, B s = B v at the same angle.
Therefore, let the incident radiation from the observation direction irradiance be 1 W/m 2 .According to the aforementioned calculation method of actual incident irradiance on the surface, if in the observation direction the irradiance on the surface is greater than zero, then the surface can be "observed".
These areas are set to 1 in the matrix B v , and the others are set to 0. Then, the binary visibility matrix B v is obtained, such as: where b rc = 0 or b rc = 1.The binary matrix B v represents the masking of the terrain in the observation direction.
Then, the reflected radiance L v in a specific observation direction can be expressed as: F. The directional-directional reflectance factor of the area observed by the remote sensing pixel Limited by the spatial resolution of remote sensing instruments, it is often difficult to obtain objects with uniform radia-tion characteristics at the pixel scale.The effect of topography on reflectance is usually concentrated on the comprehensive effect of the microslope in a remote sensing pixel.If we obtain high spatial resolution DEM data of the area corresponding to the pixel, the MRRT model can be used to obtain the directional-directional reflectance factor (BRF) (see Appendix B) of the area observed by the pixel.
Therefore, considering the effect of microtopography in a single remote-sensing observation pixel, the reflectance of microtopography is assumed to be ρ, and Lambert reflection is considered.Let T n = ρ n−1 • D n−1 , (n > 2).D n−1 represents the (n − 1)th reflection effect between terrains.When n = 1, D 0 = T 1 = E 0 .L ref can be expressed as: Then, the BRF of the observation area of a single pixel is shown in Equation (25): where θ 0 , ϕ 0 , θ v , and φ v are the solar zenith angle, solar azimuth angle, observation zenith angle, and observation azimuth angle, respectively.E sun is the direct solar irradiance.When the topography is completely horizontally flat, there is only one reflection from the terrain surface.Therefore, D 0 = E sun cos(θ 0 ), B v = 1, and BRF = ρ v .
Equation (25) establishes the relationship between the BRF of the observed pixel and the true reflectance of the microterrain inside the observed pixel area, which shows that the terrain effect on the BRF mainly comes from the true reflectance of the terrain surface, the terrain undulation, the incident irradiance on the terrain surface and the masking in the observation direction.

A. LOLA SLDEM
The Lunar Orbiter Laser Altimeter (LOLA) [34] is an instrument on the Lunar Reconnaissance Orbiter (LRO) spacecraft designed to acquire high-precision topographic data on the lunar surface [35].Barker et al. [36] combined the DEMs obtained from the LOLA DEMs and SELENE terrain camera (TC) to produce a near-global DEM with higher geodetic accuracy, namely, SLDEM2015.SLDEM2015 can be found at http://imbrium.mit.edu/DATA/SLDEM2015/.This dataset is used in this paper, and the spatial resolution of the data is 512 pixels per degree with approximately 60 m per pixel.The obvious vertical stripes obtained by visual observation of LOLA SLDEM are preprocessed by replacing the stripes themselves with the mean values of the left and right columns of the stripes.As shown in Fig. 7, the processed image is visually transitional.

B. In situ measurements from the Chang'e-3 landing site and Apollo 16 lunar soil sample 62231
The Chang'e-3 landing site is located at 44.1205 • N, 19.5102 • W [37], and a large area around the landing site has a homogeneous composition.Compared to the MS-2, Apollo 15, and Apollo 16 highland sites, the CE-3 site is much younger and less impacted and contaminated [38].Therefore, the Chang'e-3 landing site and its in situ spectra are proposed as an optical standard for both radiance calibration and wavelength calibration for lunar and Earth-orbital missions [38].
The Visible-Near Infrared Spectrometer (VNIS) instrument onboard the Yutu probe made measurements at four different points (E, S3, N203, and 205, as shown in Fig. 8) and obtained data in detection mode four times and calibration mode three times.All these data are available at https://moon.bao.ac.cn/.Apollo 16 is located at approximately 15.5 • E, 9 • S, on the relatively flat Cayley Plain, adjacent to the rugged Cartesian crater.Apollo 16 collected rock samples from the highlands, and, after that, some sample parameters were measured by the Lunar Soil Characterization Consortium (LSCC) laboratory.Therefore, the area near the Apollo 16 lunar landing site is generally selected as the research area to invert the reflectance, albedo, and mineral content of the lunar surface.Apollo 16 lunar soil sample 62231 is a typical representative of this landing site, which is characterized by high maturity and stable spectral properties [1].
The reflectance factor (REFF) spectra [38] for the CE-3 in situ measurements and the bidirectional reflectance spectra [1] of lunar soil sample 62231 are shown in Fig. 9.

C. IIM orbit data
The interference imaging spectroradiometer (IIM) onboard Chang'e-1 achieved 84% coverage of the lunar surface between 70 • S and 70 • N using push-broom hyperspectral imaging, with a spatial resolution of 200 m per pixel.The IIM achieves 32-band multispectral observations over the spectral range 480-960 nm [40].The first five bands (480-513 nm) and the last band (946 nm) of the IIM are abnormal and should be eliminated before use [38].
In this study, the data we used are CE-1 IIM 2B radiance data.The solar irradiance J at the surface of the IIM was calculated using Equation ( 26) based on the known spectral curve of the solar irradiance of the surface and the corresponding spectral response function of the IIM [17].The solar irradiance E sun is taken from the New Synthetic Gueymard Spectra [41] because it is used by the CE-3 team for the onboard calibration of the VNIS [38].
Here, f (λ, σ) is the corresponding spectral response function of the IIM, and λ 1 and λ 2 are the start and end wavelengths of f (λ, σ), respectively.The spectral response function of the IIM can be simulated with the central wavelength and full width at half maximum (FWHM) as follows [17]: In Equation ( 27), σ = F W HM 2 √ 2 ln 2 , and λ c is the center wavelength.

A. Comparative analysis of radiance for a single reflection and multiple reflections
Six regions are selected on the lunar surface (see Fig. 10).The topographic data are obtained from the LOLA SLDEM with a spatial resolution of approximately 60 m.The solar zenith angle, solar azimuth angle, observation zenith angle, observation azimuth angle, and solar incident irradiance are assumed to be 30  reflection is very abrupt, while that in the images with multiple reflections is gentle.Compared with the radiance result of a single reflection, the result of multiple reflections is more consistent with the real reflected radiance of the terrain, and more topographic details are shown.Moreover, the radiance with multiple reflections is larger than the radiance with a single reflection because the irradiance on the microtopography surface is increased by the multiple reflections between terrains.When the microtopography reflectance is small, the increase in radiance is not obvious, but as the reflectance increases, the increase in radiance is larger.This indicates that the effect of multiple reflections of radiation between terrains is nonnegligible on the lunar surface regions.
In addition, subfigures (e) and (f) in these figures all show some stripes, which are not eliminated when preprocessing the image because it is difficult to find them visually in the SLDEM image.From another view, the results show that the topography information is not broken when the calculation of the MRRT is performed, that is, the MRRT is reliable.Furthermore, the new model more easily verifies the reliability of the original data from the results than the one reflection irradiance formula, so the preprocessing of the original data can be checked again.

B. Effect of microtopography on the directional-directional reflectance of the target area
The solar zenith angle, solar azimuth angle, observation zenith angle, observation azimuth angle, and solar incident irradiance are assumed to be 30 • , 0 • , 0 • , 0 • , and 100 W/m 2 , respectively.The percentage increase in the BRF (see Equation (B2)) of the target regions relative to the microarea reflectance ρ is calculated.
Fig. 17 shows that for the same area, the difference between the BRF of the target region and the true reflectance of microtopography increases significantly with the increase in the reflectance of the microtopography.This is because the increase in microtopography reflectance will lead to an increase in multiple reflections between terrains, thus increasing the incident irradiance of the microtopography surface and the reflected radiance in the observation direction.According to Equation ( 25), the value of the denominator remains the same, and the numerator L v increases with the increase in the number of reflections between terrains, so the BRF increases.
For the same microtopography reflectance, the BRFs of different areas with the same illumination observation geometry also have obvious differences, and the differences increase significantly with the increase in microtopography reflectance.

C. Inversion of microtopography reflectance with the same illumination observation geometry
The solar zenith angle, solar azimuth angle, observation zenith angle, and observation azimuth angle are assumed to be 30 • , 0 • , 0 • , and 0 • , respectively.When the solar incident irradiances are 1, 10, and 100 W/m 2 , the reflected radiance of a completely flat region, assuming its microtopography reflectance is 0.15, can be obtained, i.e., 0.0413, 0.4135, and 4.1350 W/m 2 /sr, respectively.The calculation formula is as follows (let the flat region be Lambertian): where E f lat , ρ f lat , and L f lat are the direct solar irradiance, the microtopography reflectance of the flat region, and the reflected radiance of the flat region, respectively.θ 0 is the solar zenith angle.
Let the reflected radiance of the flat region be the observed radiance, then the microtopography reflectances of the six regions are retrieved.
Table II shows that the inversion reflectance of the microtopography of the six regions is less than 0.15 (i.e., less than the microtopography reflectance of the flat area).Solar irradiance has little effect on the inversion results of microtopography       reflectance; the order of magnitude of the inversion reflectance standard deviation for different solar incident irradiances is between 10 −5 and 10 −6 .However, under the same solar irradiance, the inversion reflectance of the microtopography of different areas has an obvious difference.
By sorting the results according to the average slope of the six areas, it can be found that the reflectance inversion results of the microtopography do not decline monotonically with increasing average slope.This is because, due to the relief of the terrain, different illumination observation geometries will cause a change in the reflected radiance in the observation direction, which is not necessarily monotonically increasing or decreasing.Therefore, if the multiple observed regions have the same illumination observation geometry and the same incident irradiance and observed radiance, the relationship between the reflectance inversion results and the topography is still unpredictable.This further demonstrates the importance of simulating multiple reflections of radiation between terrains to obtain the real terrain reflectance.

D. Inversion of microtopography reflectance by the IIM orbit data
The Chang'e-3 landing region and Apollo 16 landing region are taken as examples.We select IIM orbit data with numbers 2565, 2224, and 2843.Of these, IIM 2565 contains the Chang'e-3 landing site, and 2224 and 2843 contain the Apollo 16 landing site.The IIM orbit information is shown in Table III.
The in situ radiance and REFF (see Equation (B1)) of the Chang'e-3 landing site are taken as the standards to calculate the reflectance factor of the IIM orbit.The inversion target region uses the region shown in Fig. 10(a), which is 100×100 pixels in size and has a spatial resolution of approximately 60 m.The microarea reflectance retrieved by MRRT is the  slope surface reflectance constructed by DEM, which is the real terrain surface reflectance.Fig. 18(a) shows that the REFF of the in situ measurements is almost "5" > "6" > "7" > "8", while the REFF of calibration with the in situ measurements is "7" > "6" > "5" > "8".Moreover, the REFF of the calibration and the REFF of the in situ measurements of the corresponding site are very different.The curves of "Calibration 5" and "Calibration 8" are located among the curves of "5", "6", "7", and "8".The inversion microtopography reflectance curve by the MRRT model is located above the curve "Calibration 5" and bottom of the curve "Calibration 6", and near the curve "5".For the same region, the reflectance curves should be similar.It can be found that the curves of "Calibration 5", "Calibration 8" and "MRRT" are consistent with this inference.Among them, the inversion results of the MRRT model are larger, which may be because the reflected radiance and reflectance measured in situ at the four stations represent fewer around areas, while the Chang'e-3 landing zone selected in this paper is slightly larger than the area observed at the four stations.Therefore, there may be other microareas with large reflectance inside the target region (i.e., the Chang'e-3 landing region selected in this paper), which leads to the large reflectance inversion results of the microtopography retrieved by the MRRT model.
In addition, both calibration and MRRT curves are calculated under the illumination observation geometry of IIM radiance data itself, while the illumination observation geometry of curve "5", "6", "7", and "8" is (30 • , 0 • , 30 • ), which may also be the reason why the calibration reflectance is quite different from the in-situ measurement reflectance.However, because the MRRT model in this paper assumes that the microarea surfaces are Lambertian, the MRRT curve can still be compared with the in situ measurement curves.
Fig. 18(b) shows that the calibration REFF results of IIM 2224 and 2843 are obviously two sets, in which the set of IIM 2843 is larger than IIM 2224 at each site.Some of the calibration curves are near the curves of 62231 reflectances.The inversion reflectance curves of IIM 2224 and 2843, which have different illumination observation geometry, are very close.This proves the accuracy of MRRT reflectance inversion to a certain extent.The difference between the two curves may be because MRRT assumes that the microregion is Lambertian, while the actual terrain situation is not Lambertian in a complete sense, so there are some differences in the inversion results.Although some calibration curves are close to the 62231 reflectance curve, while the MRRT inversion results are much smaller than the 62231 curves, studies have shown a difference in composition between the Apollo 16 landing site and the actual lunar soil sample 62231, with a reflectance spectrum nearly twice as large as that observed by remotesensing sensors [14].The inversion results of the MRRT in this paper also accord with this conclusion.
For Figs. 18(a) and (b), the trend of the reflectance of all curves is similar.The last band (the 31st band of IIM) reflectance at the calibration curve and MRRT curve is sharply reduced compared with the previous band.The 31st band of the IIM may also be inaccurate.

V. CONCLUSION
In this paper, a new model for quantifying multiple reflections of radiation between terrains (MRRT) is presented.Based on the adjacent terrain irradiance formula of the first-order reflection proposed by Proy, the second-order to the nth-order reflections of radiation between terrains are derived.According to formula (25), the MRRT model establishes the relationship between the BRF of the observed pixel and the true reflectance of the microterrain inside the observed pixel area, which shows that the terrain effect on the BRF mainly comes from the true reflectance of the terrain surface, the terrain undulation, the incident irradiance on the terrain surface and the masking in the observation direction.
The establishment of the new model is helpful to the study of regional radiation characteristics of the lunar surface and to select suitable lunar surface radiation calibration fields.Compared with the ROLO model based on modeling the integration of the radiance of the entire lunar surface, the  The model presented in this paper can be used to retrieve the true reflectance of various topographic regions of the lunar surface.When using the model, it is recommended to apply it to an area with a high-resolution DEM.Generally, the higher the DEM accuracy is, the higher the accuracy of the modeling of multiple reflection processes.In addition, the MRRT model is recommended for use in combination with the Hapke model for retrieval of material in complex terrain areas; that is, the reflectance retrieved from the MRRT model is used as the true reflectance of the terrain.
Nevertheless, the model is still worth improving the accuracy of each quantity in Equations ( 16) to (25), such as the initial incident irradiance E 0 on the slope surface and the mutually visible radiation factor Γ between terrains.The research of this part can refer to the radiosity method of computer imaging, which has some similarities with the model in this paper.In addition, the application of the model to a larger scale of the lunar surface is still worthy of further study.
Our future work will focus on measuring the specific accuracy of the MRRT model in lunar calibration, such as comparison with the ROLO model, to more comprehensively evaluate the advantages of the proposed model in regional calibration.

ACKNOWLEDGMENT
We thank all the editors and reviewers for their enthusiasm and contributions to this article.They provided many valuable suggestions, which are of great help to the improvement of the content of this paper.

APPENDIX A IMAGE ROTATION ERROR UNDER DIFFERENT INTERPOLATION METHODS
In this paper, the image is rotated to match the azimuth of the solar incident radiation, so that the relative azimuth of the image and the sun is 0 • .Since the coordinates of the rotated image pixels are no longer integers, the values of the rotated pixels must be interpolated.Commonly used image interpolation methods are nearest neighbor interpolation, bilinear interpolation, and bicubic interpolation.Equations (A1) and (A2) are used to calculate the errors of the first-time rotation and second-time rotation (opposite to the direction of the first rotation) of the six groups of DEMs (see Fig. 10), respectively, and the results are shown in Table A1 and Table A2.Obviously, error 1 and error 2 are very small, among which the bicubic interpolation method has the smallest errors.
where means the sum of all elevation values in the DEM.Each rotation of the image increases the size of the image.

APPENDIX B DEFINITION OF REFLECTANCE USED IN THIS PAPER
1. Reflectance Factor (REFF): The ratio of the radiance reflected from the surface into a given direction to that of a standard panel and corrected with the REFF of the standard panel at the measurement geometry.The equation is [38]: REF F (λ, θ 0 , ϕ 0 , θ v , φ v ) = I sample (λ, θ 0 , ϕ 0 , θ v , φ v ) I std (λ, θ 0 , ϕ 0 , θ v , φ v ) × R std (λ, θ 0 , ϕ 0 , θ v , φ v ) where λ, θ 0 , ϕ 0 , θ v , and φ v are the wavelength, Sun zenith angle, Sun azimuth angle, view zenith angle, and view azimuth angle, respectively.I sample is the radiance of the target measured by the instrument, I std is the radiance of the diffuser panel measured by the instrument, and R std is the REFF of the diffuser panel.
2. Bidirectional Reflectance Factor (BRF): The ratio of the radiance reflected from the surface into a given direction to that of a perfectly diffuse surface under the same illumination observation geometry.The equation is: B2) where λ, θ 0 , ϕ 0 , θ v , and φ v are the wavelength, Sun zenith angle, Sun azimuth angle, view zenith angle, and view azimuth angle, respectively.L v is the radiance measured by the instrument, and J is the solar irradiance at the surface of the IIM, which is calculated by Equation (26).D is the Sun-Moon distance in kilometers at the observation time divided by the standard Sun-Moon distance (149,597,870 km).In this study, the value of D is 1.

Fig. 1 .
Fig. 1.Slope P reflects radiation to slope M , and P adds a source of irradiation for M .The solid yellow line is the solar incident radiation; the yellow dashed line is the absorbed solar radiation; and the solid blue line is the radiation reflected by the terrain.

Fig. 2 .
Fig. 2. Schematic representation of an "update" of a radiative source for a microarea.The dashed lines represent that there is no incident radiation, and the solid lines indicate that there is incident radiation.Arrows of different colors represent the reflected radiation from different microareas, in which the microareas indicated by arrows represent a mutually visible.In contrast, the microareas with no arrows indicate that there is no visible.The first layer receives the solar incident radiation, and the second to nth layers receive the reflected radiation from other slopes in the previous layer.

Fig. 3 .
Fig. 3.There are three types of cases in which any slope receives solar radiation: (a) fully irradiated, (b) partially irradiated, and (c) no irradiation.θ h is the solar altitude angle.

Fig. 4 .
Fig. 4. The mutual reflection between microarea surfaces is reduced to the mutual visibility between grid points.One grid point represents a flat cylinder whose height is the elevation value of the terrain.In (b), A is the radiation source, and B and C are the target receiving points; the yellow line is the radiation emitted by A; and the double circles are the actual elevation positions on the radiation path.

Fig. 5 .
Fig. 5. Schematic diagram of the stored elements in the F (1,1) matrix corresponding to the position of the image.

Fig. 6 .
Fig. 6.Schematic diagram of the actual position and elevation of the intermediate slope point.

(
iii) For max(|x a − x b |, |y a − y b |) > 1, if |x a − x b | < |y a − y b |, the number of intermediate slope points is num = |y a −y b |−1.Let the intermediate slope point be D(x d , y d , z d ).

Fig. 7 .
Fig. 7. DEM area display and image preprocessing results.(a) The original image; (b) the processed image.

Fig. 8 .
Fig.8.Map of the path traversed by the Yutu rover and the distribution of detection points[39]

Fig. 10 .
Fig. 10.Six regions on the lunar surface.The spatial resolution of the DEM in all regions is approximately 60 m.

( a )
Fig. 11.CE-3 landing site.Note that the "Proy model" refers to the result of the calculation term of the adjacent terrain irradiance in the Proy model, which is the result of a single reflection of radiation between terrains.

Fig. 17 .
Fig. 17.Increased percentage of BRF for different target areas relative to the inner microtopography reflectance with the same illumination observation geometry.The reflectance ρ values of the microtopography inside the target regions are 0.03, 0.15, and 0.3.

Fig. 18 .
Fig. 18.Inversion reflectance spectra.Legend "Calibration 5" represents the reflectance factor calibration with in-situ measurements at site 5 of the Chang'e-3 landing site, and the legend "MRRT" represents the inversion microtopography reflectance by the MRRT model in this paper.

00E+00 bilinear - 1 .
27E-02 -1.26E-02 0.00E+00 -1.27E-02 -1.26E-02 0.00E+00 -1.27E-02 -1.26E-02 0.00E+00 -1.Therefore, I, I 1 , and I ′ 2 are the original DEM, the DEM after the first-time rotation and the DEM cropped to the original DEM size after the second rotation.We align the center of the original image I and the second-time rotated image I 2 and then intercept the corresponding region according to the size of I, denoted as I ′ 2 .

TABLE II INVERSION
RESULTS OF MICROTOPOGRAPHY REFLECTANCE WITH THE SAME ILLUMINATION OBSERVATION GEOMETRY.
1The standard deviation of the retrieved reflectance for different solar incident irradiance is calculated.

TABLE III IIM
ORBITAL INFORMATION CORRESPONDING TO THE CE-3 AND APOLLO 16 LANDING SITES Area IIM orbit Instrument Zenith Angle Instrument Azimuth angle Solar Zenith Angle Solar Azimuth Angle Phase Angle 1Angles are in degrees.

TABLE A1 THE
FIRST-TIME ROTATION ERROR OF THE IMAGE: error 1

TABLE A2 THE
SECOND-TIME ROTATION ERROR OF THE IMAGE: error 2