DEM Based Study on Shielded Astronomical Solar Radiation and Possible Sunshine Duration under Terrain Inﬂuences on Mars by Using Spectral Methods

: Solar radiation may be shielded by the terrain relief before reaching the Martian surface, especially over some rugged terrains. Yet, to date, no comprehensive studies on the spatial structure of shielded astronomical solar radiation (SASR) and the possible sunshine duration (PSD) on Mars have been conducted by previous researchers. Previous studies generally ignored the inﬂuences of the terrain on the SASR and PSD, which resulted in a corresponding unexplored ﬁeld on SASR. The purpose of this paper is to study the Martian spatial-temporal structure of SASR and the PSD under terrain inﬂuences. In this paper, the theory of Earth’s SASR, the previous Martian SASR model and the theory of planetary science were combined to propose the SASR model that can be applied to Mars. Then, with the spectrum method theory of geography, we deﬁned two new concepts of spectrums to explore the spatial-temporal distribution of SASR and PSD in different Martian landforms. We found SASR and PSD on Mars were signiﬁcantly inﬂuenced by terrain relief and latitude and showed sufﬁcient regularity, which can be concluded as a gradual attenuation with terrain relief and a regularity of latitude anisotropy. The latitude anisotropy feature is a manifestation of the terrain shielding effect. With the latitude varying, SASR and PSD at different temporal scale generally showed different features with those of Earth, which may be attributed to the imbalanced seasons caused by Martian moving orbits and velocity. Compared to PSD, SASR showed more regular variation under terrain relief and was more inﬂuenced by the terrain relief which revealed that SASR is more sensitive to terrain relief than PSD. Additionally, the critical area is a quantitative index to reﬂect the stable spatial structure of SASR and PSD in different landforms and may be viewed as the minimum test region of sample areas. The corresponding result of the experiments herein indicated that either spectrum can effectively depict the spatial-temporal distribution of SASR and PSD on Mars under terrain relief and deepen the understanding of the variation of SASR and PSD inﬂuences by terrain. The critical area of either spectrum can be employed to explore and determine the stable spatial structure of SASR and PSD in different landforms. The proposed Martian SASR model and the new spectral method theory shed new light on revealing the spatial-temporal structure of SASR and PSD under terrain inﬂuences on Mars. three regions in different latitudes (the corresponding latitude area AF, CG and ED respectively). AF is the equator and ED is the Arctic Circle. When sun directly illuminates the northern hemisphere (From vernal equinox to autumnal equinox), with the increase of the latitude, the proportion of the day arc to the circumference of latitude gradually increases. Until the latitude increased to Arctic Circle, the day arc accounted for 100% of the circumference of latitude. The ratio of day arcs to latitude is the ratio of daytime to one Martian day. That is, the duration of daytime increase with the increase of the latitude of the sample region.


Introduction
Solar radiation, as an important power source, plays a vital role in agricultural production and applications of photovoltaic power generation systems on Mars [1]. photovoltaic power systems, possessing many advantages, including high power-to-weight ratio, high efficiency, modularity, scalability, and a long history of successful application in space, are required on the mission to the Martian surface [2]. Detail information on the characteristics of solar radiation on Mars is necessary for the effective design of future planned photovoltaic systems operating on the surface of Mars [3]. The exploration of the spatialtemporal distribution of solar radiation on Mars will help to deepen our understanding of effect was still out of consideration. As described, the effect of surrounding terrain may increase the huge uncertainty of solar radiation, especially on Mars, which is much more rugged than Earth. In other words, we may get an inaccurate spatial structure of solar radiation in mountains. Ultimately, due to the complex topology, the huge area, and the limitations of existing data, the existing previous research generally analyzed the solar radiation globally or at one certain site macroscopically and quantificationally, which failed to analyze the spatial distribution at a regional scale. Also, the current classification of the Martian topography is artificial. However, Mars has not been sufficiently explored. So far, only USGS and NASA have made simple and rough terrain classification, possibly resulting in an inaccurate terrain division. The specific landform has specific solar radiation characteristics, and the research on artificially defined sample area may be unable to accurately reveal the characteristics in the spatial-temporal distribution of solar radiation for the specific landform. These deficiencies lead to certain limitations of current research on SASR and PSD.
These limitations may lead us to obtain the inaccurate characteristics and spatialtemporal structure of SASR and PSD on Mars. Besides, we may get a wrong cognition of the spatial-temporal distribution variations of SASR and PSD which may be caused by the influence of the terrain. Moreover, since the terrain may influence the spatial structure of the solar energy received at the terrain surface, we may fail to be aware of how the terrain affects the received solar radiation for solar energy collectors on Mission to Mars. Given the above, under the premise of constructing the SASR and PSD model, it is of the utmost importance to find an effective method which is cable to study the SASR's spatial and temporal distribution under terrain influences as well as the influences of terrain relief effect on SASR's distributions for different Martian landforms. To this end, this research seeks to overcome these limitations.
The spectrum method, which has not been used in the study of Mars, is an effective method to reflect the spatial and temporal distribution of phenomena according to the law of gradual varies or distribution of certain indexes. It has been developed to successfully reveal the spatial structure of the slope, i.e., the slope spectrum [76][77][78][79][80]. Moreover, it was first successfully employed in SASR and PSD on Earth to reveal the spatial structure of SASR and PSD based on slope spectrum [11,12]. It is an organic combination and bold attempt of Earth science and astronomy theory. Since the Martian geomorphology meets the statistical fractal characteristic, namely, the self-similarity characteristics of the regional geomorphology required by the spectral method [79,81], the spectral method can be applied to Martian research.
The spectrum method, as an effective method, can depict the spatial structure more accurately and comprehensively than other methods [11,12,[77][78][79]. Moreover, the determination of critical areas in deriving spectra is an indispensable task [78,79]. When it is employed to analyze the SASR and PSD, it can accurately extract the stable area of SASR and PSD for different specific geomorphology [11,12], to avoid making wrong judgments on solar radiation characteristics. In other words, it is efficient in determining the stable spatial area for SASR and PSD of different Martian landforms, which may help us avoid the uncertainty caused by artificial division. Specifically, it can intuitively present the spatial distribution of SASR under different topographic relief in the form of a two-dimensional spectrum, which is convenient for us to analyze the degree of the SASR and PSD affected by topography and explore the spatio-temporal variation rules of SASR affected by topography. It is an effective method that can accurately and comprehensively reveal the distribution of SASR and PSD in specific landform and quantify the SASR and PSD under different topographic relief.
By constructing a two-dimensional spectrum of SASR and topographic relief factor, the extent of SASR affected by topographic relief can be effectively evaluated. By constructing the SASR spectrum at different temporal scales (season, month, and year), we can deepen our understanding of the variations and differences for geomorphic shielding effect of solar radiation at different durations. These advantages suggest that it can deal with most of the current existing bottlenecks mentioned above. Therefore, the spectral method is undoubtedly a valuable method for the special research status of Mars.
The specific objective of this study is to explore the spatial distribution of SASR and PSD and reveal the spatial relationship between SASR and surrounding terrain effects. The innovation points of the research can be divided into the following aspects: (1) Based on the Martian solar radiation theory, combined with Earth's SASR theory and astronomical theories, we offer the theories of SASR which have not been used on the Martian surface. This theory takes the terrain influences into account. Also, the unit for calculating the SASR is improved from one sol to twenty minutes. It somewhat enhances the shortcomings of the current Martian solar radiation model and can be employed as the basis for future Martian solar radiation calculations. (2) According to the particular topography and the limited data sources of Mars, to investigate the spatial-temporal distribution of SASR at a regional scale for different Martian landforms, one concept of spectrum method is proposed in this paper. It is an effective framework in which researchers can explore the spatialtemporal structure of SASR and PSD influenced by topography. The theory of this paper can provide a reference for the analysis on the spatial structure of other meteorological factors. Furthermore, by extracting the stable areas for different landforms, it is an efficient method that can accurately and comprehensively reveal the SASR's distribution of specific landforms regionally and quantify the solar radiation under different topographic relief. It can be concluded that the aim of this paper is to fill the gap in the studies on the analysis of spatial-temporal law of SASR on Mars, enrich and improve the theories and methods of solar radiation on the Martian surface.

Materials
The representative landforms on the Martian surface are volcanoes, impact craters, Vallis, mountain range, planitia, plateau, dune fields, and Colinas [82,83]. The analysis and interpretation of Martian geomorphologic features are mainly focused on the above mentioned typical geomorphologies [84][85][86][87][88][89]. Note that dunes are too flat to be extracted by spectral method (it may be lack of topographic feature information), thus this paper focused on the other seven landforms (further discussion see Section 4.3). All other typical landforms of Mars mentioned above which have been deeply investigated by previous researchers have been used in this paper as sample areas (see Table 1). Note that in order to confirm and demonstrate some conclusions for this paper and analyze the degree of SASR and PSD affected by shielding effect in different geographical locations, the latitude of sampling area ranges from north latitude 39.07 • to south latitude 58.1 • (see Table 1). The spectral method was successfully used in these sample areas and the stable area was extracted successfully. The data to calculate SASR and PSD is the Mars MGS MOLA -MEX HRSC Blended DEM. It is a blend of digital elevation model (DEM) data derived from the Mars Orbiter Laser Altimeter (MOLA) and the High-Resolution Stereo Camera (HRSC). It is the latest data improved from MOLA data. Its resolution accuracy is 200 m, which is significantly improved compared to the original MOLA data. It is published by USGS Astrogeology Science Center and can be downloaded from the website: https://astrogeology.usgs.gov/search/ map/Mars/Topography/HRSC_MOLA_Blend/Mars_HRSC_MOLA_BlendDEM_Global_ 200mp_v2. All the SASR and PSD of different landforms can be calculated from the MOLA data of the corresponding sample areas.

Method to Compute the PSD and SASR 2.2.1. The Law of Revolution and Rotation of Mars
Mars is a planet. In common with the Earth, Mars moves around the sun. Besides, according to Kepler's laws, it moves around sun in an elliptical orbit like the Earth. The average distance between Mars and the sun is 1.52 astronomical unit (AU). Its orbital period is 669 solar day (687 Earth days) [90]. In common with the Earth, Mars rotates constantly. A solar day(sol) on Mars, which is its rotation cycle and only slightly longer than an Earth day, is 24.6230556 h [91]. Martian eccentricity is 0.093377. Martian obliquity of rotation axis is 25.19 • [92,93]. According to the law of revolution and rotation of Mars, we can calculate the astronomical solar radiation at the top of the Martian atmosphere. Then, we can obtain the Equations of the astronomical solar radiation on the slope surface. Based on the calculation model of SASR on the sloped surface, we can reconsider the corresponding Equations and build the distributed model (or subparagraph integral mathematical model) to calculate PSD and SASR under terrain relief on Mars. Based on the existing theoretical model of solar radiation on Mars [1,2,9,10,26,27,30,32,[37][38][39] and the Astronomical theory on Mars, combining with the theoretical model of solar radiation on Earth [6,8,11,12,14,15,17,63,64,91,[94][95][96][97][98][99][100][101], the SASR model which is suitable for Mars is proposed. It can be given by the following steps:

Calculation Model of the PSD and SASR on the Sloped Surface
The SASR at the top of the atmosphere can be calculated by the following equation [2,6,27,96]: where Q 1 is the SASR at the top of the Martian atmosphere, 1 ρ is the calibration coefficient of solar-Mars distance, I 0 is the solar constant of Mars. The SASR irradiance on the sloped surface can be calculated by the following equation: where θ 1 is the solar illumination angle(or solar incident angle) on the sloped surface. cos θ 1 , which is the reduced form of Equation (14) in [1] or Equation (2.10) in paper [27], can be calculated by the following equation [1,27,63,94,96,102]: Thus, the SASR irradiance on the sloped surface can be calculated by the following equation [6,14]: In Equation (4), where δ is the declination angle; ω is the solar hour angle; U, V, W represent the characteristic parameters related to geography and terrain. Equations to compute these factors were listed as following [6,64,96]: ϕ is the latitude, α is the slope of the study point, β is the aspect of the study point. The calibration coefficient of solar-Mars distance can be given by the Equation (6) of [2] or found in [9,10,27]: where e (0.09377) is the eccentricity of Mars orbiter, L s is the areocentric solar longitude, θ is the true anomaly of Mars, and L 0 is the areocentric longitude of Mars perihelion. Note that the value of L 0 was usually evaluated at 248 • to calculate the calibration coefficient of solar-Mars distance in the previous even recent papers [2,3,10,[27][28][29][30]41,103]. However, this was assessed by a much older study [104]. Then, it was improved to a moderately accurate estimate value: 250.99 (251 • ) [55,105]. It has been employed in a few previous research works [55,56,[106][107][108][109] and employed in the Mars Climate Database (MCD), which is a database of atmospheric statistics compiled from state-of the art Global Climate Model (GCM) simulations of the Martian atmosphere. However, even now, many articles still use the 248 • as mentioned above, thereby this paper proposed it. The solar constant of Mars(or mean beam irradiance) [2]: where I solar is the solar constant of Earth(0.082 MJ · m −2 · min −1 or 1371 KW/m 2 , in this paper, we used 0.082 MJ · m −2 · min −1 ) [2,6], R Mars is the average distance from Mars to the sun; R solar is the average distance from Earth to Sun [2]. The declination angle can be calculated by the spherical triangle law of cosines: where ε (25.19 • ) denotes the obliquity to orbit. According to [2], it can be simplified even further to: However, since the rise and set of the sun, each Martian sol has a limited amount of PSD. and We present the PSD as an interval [ω r , ω s ].. Adding up the SASR for the PSD from sunrise to sunset, we can obtain the Equation for calculating one Martian sol's quantity of SASR on the sloped surface which is the integration of Equation (4) [6,14]: where ω r1 is the solar hour angle at the begining of the PSD, ω s1 is the solar hour angle at the end of the solar hour angle, Q 0 denotes the daily SASR on the sloped surface from ω r1 to ω s1 ; T is the time of one Martian sol, i.e., 24.6230556 × 60(in minutes). ω s1 and ω r1 , i.e., the true sunrise and sunset on the incline plane can be obtained from equation 13-21 in [1]. The true beginning of the PSD and the true hour angle at the end of the PSD take into account of the hour angle of sunrise and sunset on the horizontal plane and the hour angle of sunrise and sunset on the inclined plane. The true sunrise should be determined by which of the two sunrises occurs later and the actual sunset should be determined by which of the two sunsets occurs earlier [1,27,32].
However, this paper did not use the method of [1] to determine the true sunrise and sunset of PSD. The actual surface in DEM is rugged, not stationary incline (there is terrain relief in all directions). The rugged terrain may lead to the shielding effect. As we know, both the sun elevation angle and the sun azimuth angle vary with the time (solar hour angle), which cause the sun angle of incidence is different. Therefore, in each solar hour angle, SASR may be shielded by the terrain relief. The distributed model (or called subparagraph integral mathematical model) was generally adopted to calculate the actual SASR and PSD over the rugged terrain [6,8,11,12,14,15,17,64,96,97,109] in DEM. It was given in the following part.

Theoretical Equation for Calculating PSD and SASR on Rugged Terrain
SASRs are affected by the terrain relief at the total period of PSD. That is to say, in some subperiods of PSD for certain Martian sol, the study point may be shielded by the terrain relief and the subperiods of PSD should be zero. Thus, we should reconsider Equation (13). the general method in digital terrain analysis [6,8,11,12,14,15,17,64,96,97] is that we take horizontal sunrise and sunset angle as a value domain. Then, we cut this range into n small intercell. The 20 min were used as an intercell. We determine the shielding effect of the given position each 20 min. There are two methods to determine the shading status of the sample point in each intercell. One is by the ray tracing method [8] and the other is by the hillshade algorithm [109]. Then, we can get m periods of valid PSD: [ω sr1 , ω ss1 , . . . . . . , ω srl , ω ssl , . . . . . . , ω srm , ω srm ](l = 1, 2, . . . m). The true sunset and sunrise should be the beginning of the first available PSD intercell and the end of the last valid PSD intercell. The general steps are given in the following: (1) ω r and ω s are the hour angle of sunrise and the hour angle of sunset on the horizontal plane, respectively. Their values, changing with the varies of points' location and time, can be given by: Note that the above Equation applies only to |ϕ| <π/2−|δ|. There are two special cases in which the above Equation is not applicable. Correspondingly, when − tan ϕ tan δ> +1, The phenomenon of the polar night will occur, i.e., PSD should be zero. likewise, solar radiation should also be zero; When − tan ϕ tan δ < −1, the phenomenon of the polar day will occur, i.e., PSD should be all day [1,2]. (2) Using ∆ω as a sun hour angle step length, the [ω r , ω s ] should evenly be divided into n intercell. For each intercell, we consider whether solar radiation is covered by surrounding terrain at the study point. Besides, there is a time step length denoted by ∆t for each different ∆ω. ∆ω can be given by: n can be given by: where int ( ) is a function which can obtain the integral part of the ω s −ω r ∆ω Thus n intercell can be obtained: (3) Whether a point can be exposed to sunlight depends on the sun elevation angle and the azimuth angle and the shielding status function caused by surrounding terrain.
The sun elevation angle and azimuth angle of every intercell can be given by the following Equations: where h i is the sun elevation angle; A i is the azimuth angle; Hillshade is a function of the terrain shading status depending on which we can get the terrain shading factor for each different period. we denote terrain shading status factor as S i where S i = 0 indicates the study point was shaded and S i = 1 indicates It was exposed to the sun. The hillshade algorithm, which was derived from and employed in the ArcGIS tool, was used to analyzes the local terrain shadow by considering the illumination source angle (solar incident angle) and surrounding terrain shadows. It is generally used to determine the assumed brightness value of a given position at a given time (solar altitude angle and solar azimuth angle vary with solar hour angle) influenced by the surrounding terrains. It has been used to determine the shielding status of solar radiation. [109][110][111][112][113][114][115][116][117]. The obtained value, which was named as hillshade value, ranges from 0-255. The corresponding Equation was The hillshade value ranges from 0-255. When hillshade>0, the place is incident by sun; when hillshade=0, the place is shielded by the terrain. In this paper, the hillshade was reclassified to terrain shading status factor as S i according to Equation (22). When hillshade > 0, S i = 1; when hillshade = 0, S i = 0. (4) We now can calculate the PSD with considering the shielding effect. The PSD considering shielding effect should be obtained by adding up all the PSD subperiod which has available sunshine. Defining g i as the shielding coefficient of each period which can be determined by the surrounding S i : The PSD of certain study point in the rugged terrains can be given by: where mod () is a function which can obtain the remainder part of the ω s −ω r ∆ω . Hereafter, the PSD indicates the PSD which considers the terrain relief. (5) The SASR which considers the terrain relief effect is obtained by adding up the SASR in all available PSD. Whether a point can be exposed to sunlight in the integral part [ω i−1 , ω i ] is determined by the S i−1 and S i . If S i = 0, S i+1 = 0 or S i = 0, S i+1 = 0, nothing should be done in the [ω i , ω i+1 ]; If S i = 0, S i+1 = 1, the angle of new PSD should begin as ω srl which was the mean value of ω i and ω i+1 ; If S i = 1, S i+1 = 0, the angle of new PSD should end as ω ssl which was the mean value of ω i and ω i+1 . Thus, an array of solar hour angles for an m segment can be obtained (in this array, the true sunrise and sunset solar hour angle over the rugged terrains should be the beginning of the first valid PSD and the end of the last valid PSD respectively): [ω sr1 , ω ss1 , . . . . . . , ω srl , ω ssl , . . . . . . , ω srm , ω srm ](l = 1, 2, . . . m) We can obtain the daily SASR at the calculated point by adding up the SASR in the PSD which has available sunshine: where W 0 denotes the SASR for one certain Martian sol, other parameters see above.

Calculation of PSD and SASR for Each Martian Sol, Season, and One Year
To calculate SASR on Mars for different Martian days, we should compute the declination δ and areocentric longitude L s for each Martian sol. Two methods to obtain the relationship between the Martian sols and the areocentric solar longitude were given in this paper. They were given as follows: Method 1: the first method was based on Kepler's equation. It was firstly proposed by Johannes Kepler in 1609 to reflects the function relationship between the position of the planet on the planetary orbit and the corresponding solar time. For a given sol number T n , the corresponding value of areocentric solar longitude L s maybe calculated by Kepler's equation. To recapitulate briefly, it can be calculated by the following two steps: Step 1: Establishing the Kepler' s equation as: where M is the mean anomaly, E is the eccentric anomaly. The mean anomaly M is related to sol numbers by: where T i is the sol number, T s is the number of sols in a Martian year (668.6), T p is the perihelion date: 485.35. Then, under the Equation (28) and (29), for each T i , the corresponding E can be obtained. Noted that the solutions were usually dissolved by Newton iterative method which may yield an identical solution to Kepler's solution.
Step 2: calculations of the true anomaly: Geometrically, the eccentric anomaly E and the true anomaly θ have the following relationship: Then, the relationship between θ and E can be further simplified to: Step 3: calculations of the areocentric solar longitude L s : Method 2: the second method was based on Kepler's second law: a line joining a planet and the Sun sweeps out equal areas during equal intervals of time. It can be denoted by the Equation below: where a is the Mars semi-major axis, b is the Mars semi-minor axis. Other parameters see the papers above. Then, for each T i , the corresponding θ can be obtained by the above Equation (34). Under the Equations (33) and (34), the corresponding L s can both be calculated. For the above two methods, we calculated separately, and the result seemed little different (the ratio of result for the first method to the second method was between 0.999871-1.00081). In this paper, we adopted the first method to calculate the corresponding L s for each Martian sol.
The duration division on Mars, shown in Table 2, was different from the Earth due to the difference of orbit and speed between Mars and Earth. In this paper, referring to the MCD' s standard, the areocentric solar longitude was equal divided into 12 classes. Each areocentric solar longitude class corresponds to a Martian month and every three months corresponds to a season as shown in Table 2. For each sample area, using twenty minutes as a unit, we calculated the value of PSD and SASR for each twenty minutes in one Martian sol. Then, we obtained the PSD of each Martian sol by adding up the SASR and PSD of all the corresponding twenty minutes in one sol. Similarly, the PSD and SASR in each season were obtained by adding up the PSD and SASR of the days in each season. Ultimately, we can obtain the SASR and PSD for one year by adding up the SASR and PSD for 669 Martian sols. Slope has been proven to be a parameter which is a significant methodology in revealing features of SASR on Earth to some extent [12]. It is also an efficient method to reveal the relations between terrain and SASR [11,24]. However, since the limitations of MOLA data (the latest available data is not a high-precision one), the microcosmic terrain factor (i.e., slope) may be not an available factor to be employed on Mars for revealing the influences of terrain relief to SASR in different topography. A low resolution of the data will result in lower mean value and standard deviation of slope information [118]. That is, if we use slope to describe terrain relief in the resolution of 200 m, the terrain feature information may be single and inaccurate. This was also verified by our experiments. We selected 50 random sample areas to calculate their slopes. For each sample areas, the extracted average slope is generally less than 10 degrees and standard deviation is generally less than 8, which proves that slope is not suitable for the study of the terrain analysis of SASR on Mars. Accordingly, we should find an efficient macro terrain factor to replace the slope to reveal the terrain relief of Mars. In topographic statistical analysis, roughness was generally used together with slope to describe terrain relief. It may be a proper parameter to displace the slope which has been employed in many planetary theories to reveal the terrain relief including Mars and the moon [73,81,[119][120][121][122][123][124][125].
There are numerous representative models to quantize the roughness on the surface. They can be divided into surface area ratio model [126], vector ruggedness model [127][128][129], surface roughness factor based on the delta transformation of slope and aspect [129,130], commonly reported models based on statistical theory (the Hurst exponent, root-meansquare height, root-mean-square slope, height distribution deviation, autocorrelation length, median and absolute slope, etc.) respectively [125,130]. The method to quantitatively characterize surface roughness will differ according to different application purposes. The quantification of local topographic relief is usually simulated by using statistical theoretical parameters. In the statistical analysis of roughness, root-mean-square height (RMS) was a high-frequency parameter. RMS, which was generally considered as a proper parameter to describe roughness in statistical analysis, were used to produce the global roughness surface maps of Mars and analyze local terrain relief on Mars and the moon [81,119,121,[123][124][125]. It seemed to show a good ability to reveal the local terrain relief on the surface [118][119][120]123,131]. It was also used to classify different geological units [73,122]. Furthermore, it has proven to have a good ability to reveal characteristics of slope variation and local terrain variation [132]. Thus, in this paper, we use RMS to quantize the local terrain relief. It can be calculated by the following Equation: where n is the number of sample points in the local analysis window, h(x i ) is the height of the center sample point of the analysis window, and h mean is the average height of all the observed sample points. In general, the rugged the surface, the larger RMS is.
Since the RMS has scale dependence, the choice of the right step size is critical. Note that the appropriate scale varies according to our algorithm for SASR affected by topography. The terrain algorithm can be divided into in hillshade algorithm which is both used in Arc/Info and Genasys [109][110][111][112][113][114][115][116] and the ray-tracing method. When using the ray-tracing method, the appropriate scale of the calculated terrain factor should vary with the search radius of the ray-tracing method. In this paper, considering that the data resolution is 200 m (in meter), which is not the high-resolution data, we adopted the hillshade algorithm. On the other hand, the hillshade algorithm determines the shielding effect of each grids according to the eight points around the grid points [109,113,116]. That is, when using hillshade algorithm, we consider the terrain relief of the 8 sample points around the sample points, so the local terrain relief revealed by our topographic factor should also be the terrain relief around the sample points exposed. Thus, when we use roughness to reveal the local terrain relief, the proper scale for roughness should be consistent with the analysis sample window of the hillshade algorithm. The analysis window of hillshade algorithm is a 3 × 3 grid (each grid is a 200 × 200 m grid). i.e., the proper scale for roughness is the range of 600 m × 600 m.

Extracting Procedure of the Stable Roughness-Mean Shielded Astronomical Solar Radiation Spectrum
The SASR spectrum can be viewed as a group of SASR histogram under certain grading laws based on one certain terrain factor. In this paper, we used RMS as a parameter in the x-axis; i.e., the SASR spectrum was a group of SASR histogram under certain grading laws based on RMS. Likewise, the PSD spectrum is a group of PSD histogram under certain grading laws based on RMS. A specific geomorphic type has its own characteristic SASR and PSD spectrum curves. The validity of the spectrums lies on whether the spectral sample region of the spectrum we used to extract is greater than the critical area threshold. When the spectral sample region of the sampling area is greater than the critical area threshold, the characteristic spectral curve shows a stable trend, and vice versa (see Figures 1 and 2). The spectrum whose spectral region is greater than the critical area is considered as the stable roughness-mean SASR spectrum of the landform and the stable roughness-mean PSD spectrum of the landform, respectively. For convenience, the stable roughness-mean SASR spectrum may be called RS hereafter. Furthermore, the critical area of RS may be called X 1 . Similarly, the stable roughness-mean PSD spectrum may be called RP in the following paper for brevity. Then, the critical area of RP may be called X 2 . The technical route to extract RS and RP is given in Figure 3. Similarly, we can extract the stable roughness-mean PSD spectrum and obtain the critical area of (i.e., 2 ).

Figure 1.
The gradual expanding process of the spectrum. Grid A is the center of the analysis window. This yellow window consist of 5 × 5 grids is the initial analysis window. Then, it expends the periphery area to be a new rectangular window whose shape is 7 × 7. When the quantitative indicators of similarity do not meet the abovementioned conditions, it will continue to expand. As we can see, it is a gradual process: The window size goes from 5 × 5 to 7 × 7 to 9 × 9…. Figure 2. 80 random grid points were selected to yield the spectrums in Valles Marineris. When the area of the extracted analysis window was less than the critical area (the green sample window), the spectrums showed the characteristics of the unstable disorder. Afterwards, along with the increase of analysis windows size, the spectrums gradually show a certain similarity and stability. Finally, when the area of the extracted analysis window was greater than the threshold of the critical area (the black sample window), the spectrums all showed a stable trend. Besides, the extracted analysis windows whose area is greater than the threshold of the critical area were all regarded as the stable area, which is because the extracted spectrums generally showed a certain similarity and stability.

Figure 1.
The gradual expanding process of the spectrum. Grid A is the center of the analysis window. This yellow window consist of 5 × 5 grids is the initial analysis window. Then, it expends the periphery area to be a new rectangular window whose shape is 7 × 7. When the quantitative indicators of similarity do not meet the abovementioned conditions, it will continue to expand. As we can see, it is a gradual process: The window size goes from 5 × 5 to 7 × 7 to 9 × 9 . . . . Similarly, we can extract the stable roughness-mean PSD spectrum and obtain the critical area of (i.e., 2 ).

Figure 1.
The gradual expanding process of the spectrum. Grid A is the center of the analysis window. This yellow window consist of 5 × 5 grids is the initial analysis window. Then, it expends the periphery area to be a new rectangular window whose shape is 7 × 7. When the quantitative indicators of similarity do not meet the abovementioned conditions, it will continue to expand. As we can see, it is a gradual process: The window size goes from 5 × 5 to 7 × 7 to 9 × 9…. Figure 2. 80 random grid points were selected to yield the spectrums in Valles Marineris. When the area of the extracted analysis window was less than the critical area (the green sample window), the spectrums showed the characteristics of the unstable disorder. Afterwards, along with the increase of analysis windows size, the spectrums gradually show a certain similarity and stability. Finally, when the area of the extracted analysis window was greater than the threshold of the critical area (the black sample window), the spectrums all showed a stable trend. Besides, the extracted analysis windows whose area is greater than the threshold of the critical area were all regarded as the stable area, which is because the extracted spectrums generally showed a certain similarity and stability. When the area of the extracted analysis window was less than the critical area (the green sample window), the spectrums showed the characteristics of the unstable disorder. Afterwards, along with the increase of analysis windows size, the spectrums gradually show a certain similarity and stability. Finally, when the area of the extracted analysis window was greater than the threshold of the critical area (the black sample window), the spectrums all showed a stable trend. Besides, the extracted analysis windows whose area is greater than the threshold of the critical area were all regarded as the stable area, which is because the extracted spectrums generally showed a certain similarity and stability.
ISPRS Int. J. Geo-Inf. 2020, 9, x FOR PEER REVIEW 15 of 42 Figure 3. Technical route to extract the roughness-mean shielded astronomical solar radiation spectrum and roughnessmean possible sunshine duration spectrum. The roughness was reclassified by root-mean-square (RMS). The declination and solar constant for each day order calculated from Kepler's law can be employed to calculate the calibration coefficient of solar-Mars distance (see Equation (8) and Equation (9)) which is the intermediate value of the calculations for shielded astronomical solar radiation (SASR) and possible sunshine duration (PSD) on rugged terrains. Then, the final Equation (see Equation (27) More tests, presented in Figure 2, were conducted to verify the feasibility of the stable . Technical route to extract the roughness-mean shielded astronomical solar radiation spectrum and roughnessmean possible sunshine duration spectrum. The roughness was reclassified by root-mean-square (RMS). The declination and solar constant for each day order calculated from Kepler's law can be employed to calculate the calibration coefficient of solar-Mars distance (see Equations (8) and (9) That is, the SASRs of the grids whose RMS was in the same class were averaged to be one value. Then, the corresponding 17 values obtaining from each class became a group. The initial spectrum can be denoted as a group of histogram: RS 1 = {S 1 , S 2 , S 3 , S i , . . . , S n }, where i is the number of roughness classes and S n denoted the mean SASR values of grids belonging to the same RMS class.
The spectrum method reveals the similarity of nearby terrains. We defined the quantitative indicators of similarity as & 1 where RS 1 denoted the spectrum before expansion, RS 2 denoted the spectrum after expansion. & 3 is an auxiliary parameter, which requires that the grids in the spectral sample region have a least one value in each RMS class. They were employed to determine whether the roughness-mean SASR spectrum before expending is stable.
The process of obtaining the critical area is an expanding process. The gradual expanding process of the spectral sample region is presented in Figure 1. The initial spectral sample region was a rectangular analysis window of 5 × 5. We chose a random grid in the sample area as the center of the analysis window. The roughness parameter (RMS) of each grid in the analysis window was calculated as mentioned in Section 2.3.1 and classified into 17 classes. The SASR of each grid in the spectral sample region was also calculated as the method mentioned in Section 2.2. Then, we added up the SASR values of grids belonging to the same RMS class and calculate the mean value, respectively. The initial spectrum, i.e., a group of the histogram in a coordinate whose x-axis is the classified RMS and y-axis is the mean SASR of grids belonging to the same RMS class, was obtained.
The initial spectral sample region of specific landform was constantly expanded, and the quantitative indicators were constantly calculated until there are thirty continuous expanding spectrums that satisfy the & 1 < 0.001, & 2 < 0.001, & 3 . When the continuous thirty expanding spectrums satisfied the corresponding conditions, we regarded that the spatial structure of spectrums from the first spectrum tends to be stable (i.e., under the limitations of & 1 and & 2 , the corresponding spectrums derived from adjacent locations reach the similarity of 99.9%) The spectrum which was constructed by the first spectrum was regarded as RS. The area of the current spectral sample region based on the first spectrum was considered as the critical area of RS (i.e., X 1 ).
However, if it does not meet the conditions, the current spectral sample region should continue to gradually extend until & 1 , & 2 and & 3 meet the conditions. That is, when the area of the spectral sample region increased, the spectrum constructed by the current spectral sample region will be closer to the stable spectrum.
After the process of expanding under the limitations, we can obtain the RS, which is a group of the histogram: RS = {S 1 , S 2 , S 3 , S i , . . . , S n }, where i is the number of roughness classes and S n denoted the mean SASR values of grids belonging to the same RMS class (as seen in Figure 2).
Similarly, we can extract the stable roughness-mean PSD spectrum and obtain the critical area of RP (i.e.,X 2 ).

Extracting Procedure of the Standard Stable Roughness-mean Shielded Astronomical Solar Radiation Spectrum
More tests, presented in Figure 2, were conducted to verify the feasibility of the stable roughness-mean SASR and PSD spectrum. Referring to the method of [77,79], random functions were taken to select 80 grids in all sample areas as the center of the analysis windows. These grids were used to yield 80 roughness-mean SASR spectrum and roughness-mean PSD spectrum, respectively. As presented in Figure 4, these spectrums showed the same trend and stability, which prove the stability and self-similarity of the extracted stable spectrums to some extent. In this paper, the coefficient of variation (COV) and average relative error (ARE) were used to reflect the statistical dispersion of 17 RMS classes. Furthermore, the correlation coefficient (CORREL) and covariance were used to reflect the similarity of spectral curves. The calculation of correlation coefficient and covariance were not discussed in this paper due to the limited space. We counted 17 RMS class values of 80 spectral curves respectively. The calculation Equation of the COV and ARE is as follows: where j is the number of RMS class, i is the number of spectrum curve. The further the coefficient of variation is from zero, the greater the statistical dispersion of the data is.  Noted that the grid with RMS which were greater than 80 m account for less than 6% of the total area in all sample areas, so we treat the interval which is greater than 80 m as a separate class. The stable spectrums, extracted in 80 regions separately, were presented in the figure. All spectral curves are distributed in a band and have the same trend and distribution. This showed a self-similarity and stability of the spectrums. Due to the paper's space constraints, only the spectrums of Valles Marineris were given in the figure.

Characteristics of Annual Spectrums
The  Given the local spectrum may not reflect the overall trend of SASR distribution in the sample area, to assess the extent of SASR and PSD affected by the terrain relief of the whole test area in a statistical sense, the standard SASR and PSD spectrum should be considered [76][77][78][79][80]. The standard spectrum was generally defined as the average of the extracted spectrum. In this paper, we defined the standard roughness-mean SASR spectrum as R S S. R S S was obtained by averaging the extracted 80 RS. The RS in each sample area was denoted as RS = S 1,j , S 2,j , . . . S i,j . . . , S i,j where i is the number of RMS classes, j is the number of divided areas. As such R S S can be expressed as the following equation: Similarly, we defined the roughness-mean PSD spectrum as R S P and calculated it by averaging the extracted 80 RP.
We defined the standard X 1 as X 1s and the standard X 2 as X 2s . They were the maximum extracted critical area of the extracted 80 RS and RP respectively.

Characteristics of Annual Spectrums
The extracted R S S and R S P for 12 sample areas in a year are shown in Figures 5 and 6, and the corresponding X 1s and X 2s are presented in Figure 7. R S S and R S P showed the specific spatial structure of SASR and PSD for each sample area. They were generally affected by the surrounding terrain and the position of latitude, which can be divided into the following two respect to be discussed.    The blue curve is the former. The red curve is the latter. As seen in the figure, the critical area of R S S was generally greater than that of R S P.

The Law of Gradual Variation of R S S and R S P with Terrain Relief in Annual Spectrum
Result 1. We found that, for the same landforms, within each RMS class in R S S, the value of SASR was generally different, which suggested that SASR showed different characteristics in different topographic conditions under the same landform. This law was also found in the PSD in R S P.
Besides, in R S S and R S P, along with the increase of RMS classes, the value of SASR and PSD showed the same trend, i.e., they were generally decreased (see Figures 5 and 6). This can be explained by the following. Roughness is a parameter that can reveal the terrain relief of the given area. The greater the roughness, the more rugged the terrain. Then, the more rugged terrain may lead to a greater shielding effect than the less rugged terrain especially the flat surface. The shielding effect may generally cause the PSD and SASR to decrease. This finding seemed to be similar to the features of SASR and PSD found on Earth [6,8,12,14,64]. And it provided deep support for the conclusion: the region possessing complex terrain may generally have a strong terrain shielding effect.
That is, the SASR and PSD were influenced by terrain relief, either on Mars or Earth.

The Law of Critical Area with Terrain Relief
Result 2. We investigated the value of X 1s and X 2s found in 12 sample areas, which seemed to show the same descending order (see Figure 7). Furthermore, the value of X 1s was generally greater than X 2s (see Figure 7). That is, it seemed that a certain of similarities exist on the occurrence of the S 1 P and S 2 P; and the variation of SASR is larger than PSD so that a larger critical area is needed to ensure the stability of the obtained spectrum.
Whether the critical area can be obtained depends on whether the landform has enough topographic feature information. For complex landforms, the topographic feature information will also be more complex, which further cause the corresponding spatial distribution of SASR to become more complex. The more complex spatial distribution of SASR may bring more topographic feature information, which leads to the larger X 1s . Accordingly, we may draw an inference: the value of critical area may be related to the complexity of topography. The variation coefficient of RMS described the variation of the terrain relief. Then, we found that a positive relationship seemed to exist between the value of X 1s and the variation coefficient of RMS (see Figure 8). This feature was also found in X 2s (see Figure 9). Furthermore, the following Equations were given in the test of 12 sample areas: where x is the value of the variable coefficient of the RMS in each sample area and R 2 the correlation coefficient. The above Equations may reveal that the value of the critical area is generally positively correlated with the variation coefficient of RMS. The variation coefficient of RMS described the variation extent of the terrain relief. When it increased, the terrain relief may become complex. Then, the complex terrain relief may lead to a complex shielding effect which may cause the distribution of SASR and PSD more complex. Then, a larger critical area is needed to ensure the stability of the obtained spectrums. Namely, the more complex terrain generally meant the higher X 1s and X 2s .
3.1.3. Latitude Anisotropy Characteristics Influenced by the Shielding Effect in R S S and R S P Result 3. Afterwards, we investigated that the R S S usually showed an intimate relationship with the latitude of the site. As seen in Figure 5, for R S S in the sample areas, with the latitude of the sample area decreasing, the SASR for each RMS class increased. The closer the sample area to the equator, the greater the SASR for each RMS class. This characteristic can be explained by the following. The lower the latitude, the higher the altitude angle of the annual mean midday sun (explained in Figure 8), and thus the weaker the shielding effect of the sample area. The weaker shielding effect may result in greater SASR. hemisphere, the PSD showed the contrary characteristics. That is, for in the sample area of the northern hemisphere, with latitude increasing, the PSD generally increased. This characteristic may be because of the imbalanced seasons caused by Martian moving orbits and velocity (see Table 1). The unique Martian orbits and velocity influenced its revolution and rotations, which result in the corresponding time laws. It may be explained as follows. Because of the obliquity of the ecliptic, the sun can illuminate only half of the Earth, which was divided by the terminator. When the sample area in the northern hemisphere is in spring and summer, the daytime which is from sunrise to sunset may become longer with the latitude increasing (see the corresponding explanation in Figure. 9) Besides, with latitude increasing, the shielding effect may gradually enhance. Generalizing these two factors, with latitude increasing, the PSD in spring and summer may increase. Similarly, we may infer that with latitude increasing, the PSD in the northern hemisphere in autumn and winter may become decrease. Then, the unique Martian movement causes the days of spring and summer on Mars are much longer than that of autumn and winter (the ratio of the durations for each season is 1.34: 1.25:1: 1.08, while it is nearly 1:1:1:1 on Earth), which may lead to the result that the increase of PSD in spring and summer is greater than the decrease of PSD in autumn and winter. Thus, the PSD in the northern Martian hemisphere in a year, which was calculated by adding up the PSD of four seasons, may increase with the increase of latitude. Figure 9. Variations of the duration of daytime with the increase of latitude in the northern hemisphere. AO, CB and DE respectively represent the day arcs at three regions in different latitudes (the corresponding latitude area AF, CG and ED respectively). AF is the equator and ED is the Arctic Circle. When sun directly illuminates the northern hemisphere (From vernal equinox to autumnal equinox), with the increase of the latitude, the proportion of the day arc to the circumference of latitude gradually increases. Until the latitude increased to Arctic Circle, the day arc accounted for 100% of the circumference of latitude. The ratio of day arcs to latitude is the ratio of daytime to one Martian day. That is, the duration of daytime increase with the increase of the latitude of the sample region.

Figure 8
Additionally, the above result seemed to show that SASR is more sensitive to the shielding effect than PSD because it showed sufficient regularity under the influence of the shielding effect. Figure 9. Variations of the duration of daytime with the increase of latitude in the northern hemisphere. AO, CB and DE respectively represent the day arcs at three regions in different latitudes (the corresponding latitude area AF, CG and ED respectively). AF is the equator and ED is the Arctic Circle. When sun directly illuminates the northern hemisphere (From vernal equinox to autumnal equinox), with the increase of the latitude, the proportion of the day arc to the circumference of latitude gradually increases. Until the latitude increased to Arctic Circle, the day arc accounted for 100% of the circumference of latitude. The ratio of day arcs to latitude is the ratio of daytime to one Martian day. That is, the duration of daytime increase with the increase of the latitude of the sample region.
As presented in Figure 6, the same characteristics were also shown in R S P of the sample area in the southern hemisphere. However, we found that in the northern hemisphere, the PSD showed the contrary characteristics. That is, for R S P in the sample area of the northern hemisphere, with latitude increasing, the PSD generally increased. This characteristic may be because of the imbalanced seasons caused by Martian moving orbits and velocity (see Table 1). The unique Martian orbits and velocity influenced its revolution and rotations, which result in the corresponding time laws. It may be explained as follows. Because of the obliquity of the ecliptic, the sun can illuminate only half of the Earth, which was divided by the terminator. When the sample area in the northern hemisphere is in spring and summer, the daytime which is from sunrise to sunset may become longer with the latitude increasing (see the corresponding explanation in Figure 9) Besides, with latitude increasing, the shielding effect may gradually enhance. Generalizing these two factors, with latitude increasing, the PSD in spring and summer may increase. Similarly, we may infer that with latitude increasing, the PSD in the northern hemisphere in autumn and winter may become decrease. Then, the unique Martian movement causes the days of spring and summer on Mars are much longer than that of autumn and winter (the ratio of the durations for each season is 1.34: 1.25:1: 1.08, while it is nearly 1:1:1:1 on Earth), which may lead to the result that the increase of PSD in spring and summer is greater than the decrease of PSD in autumn and winter. Thus, the PSD in the northern Martian hemisphere in a year, which was calculated by adding up the PSD of four seasons, may increase with the increase of latitude.
Additionally, the above result seemed to show that SASR is more sensitive to the shielding effect than PSD because it showed sufficient regularity under the influence of the shielding effect.
Result 4. Another investigation is that the descending rates of SASR in R S S and PSD in R S P are affected by latitude (see Figures 10 and 11). In R S S and R S P, the relationship between SASR, PSD and roughness classes may be generally denoted as a linear relationship: y = a n x + b, which y denotes the value of SASR or PSD, x denotes the number of roughness classes and the R 2 is ranged from 0.9012-0.998. a n is generally negative, representing a decreasing trend of SASR and PSD with the increase of roughness classes, which may provide support for result 1. Then, we found that, as seen in Figures 10 and 11, along with latitude decreasing, the a 1 and a 2 for SASR and PSD respectively, seem to generally increase. The closer the sample area to the equator, the greater a n . The explanation for this conclusion is an inference from the previous conclusion. According to result 2, with the decrease of latitude, the sun elevation angle may increase, resulting in a weaker shielding effect in S s P. Besides, the decay rate of the SASR (a 1 ) and PSD curve (a 2 ) represented a variation rule of SASR and PSD under the influence of topographic relief. As the shielding effect weaken, the PSD within each RMS class will approach the ideal value of PSD (the latter did not consider the shielding effect and was a constant value, i.e., the PSD in the length of the daytime from sunrise to sunset). The smaller the latitude, the closer the values in each RMS class are to the same value. That is, as the latitude decreases, a 2 may become gradually greater, which represents the continuous decrease of the decay rate. The explanation for the increase of a 2 with the decrease of latitude can be given in a similar way provided for a 1 .
Furthermore, we found that a 1 was generally less than a 2 (see Figures 8 and 9). It meant that the attenuation rate of PSD with topographic relief is lower than that of SASR under the same terrain shielding condition. This finding may further provide support for the inference in result 3 that SASR is more sensitive to the shielding effect than PSD.
crease of the decay rate. The explanation for the increase of 2 with the decrease of latitude can be given in a similar way provided for 1 . Furthermore, we found that 1 was generally less than 2 (see Figure 8 and 9). It meant that the attenuation rate of PSD with topographic relief is lower than that of SASR under the same terrain shielding condition. This finding may further provide support for the inference in result 3 that SASR is more sensitive to the shielding effect than PSD. The descending rate of SASR and PSD showed the same trend: from the northern hemisphere to the equator，the rate became slower; Then, from the equator to the southern hemisphere, the rate becomes faster.

Characteristics of Spectrums in Four Seasons
Since Martian moving orbits and velocity are different from Earth, the durations of Martian seasons are not evenly distributed. This paper divided four seasons of Mars according to table 2(employed in MCD). Then, we extracted and in four seasons respectively. The employed method was in common with the method used in the annual spectrum.  The descending rate of SASR and PSD showed the same trend: from the northern hemisphere to the equator, the rate became slower; Then, from the equator to the southern hemisphere, the rate becomes faster.

Characteristics of Spectrums in Four Seasons
Since Martian moving orbits and velocity are different from Earth, the durations of Martian seasons are not evenly distributed. This paper divided four seasons of Mars according to Table 2 (employed in MCD). Then, we extracted S s R and S S P in four seasons respectively. The employed method was in common with the method used in the annual spectrum.
As seen in Figures 12 and 13, since the spectrums in four seasons represent stable spatial structures of SASR or PSD for four different seasons on Mars, the seasonal combination spectrum can be used to reveal the temporal distribution of spatial structure for SASR and PSD in different landforms under terrain relief. We used it to analyze the variation law of SASR and PSD in different sample areas under topographic relief.

The Law of Gradual Variation of R S S and R S P with Terrain Relief in the Annual Spectrum
The spectral curve of S s R and S S P in each season generally showed a downward trend with the increase of the number of RMS classes in each season, which was similar to the annual spectrums (see Figures 12 and 13). The explanation for these characteristics may be consistent with the above explanation provided for the annual spectrums. That is, no matter the temporal scale, year or season, SASR and PSD keep a regular downward trend with the increase of terrain relief. We may infer that the rugged terrain generally obtained the weaker PSD and SASR due to the shielding effect caused by the more complex terrain.

Latitude Anisotropy Characteristics Influenced by the Shielding Effect in R S S and R S P
Furthermore, we found that with the latitude varying, characteristics discovered in the S s R and S S P between four seasons showed a certain of regularity. This finding will be of interest that the variation of SASR and PSD on Mars was different from that on Earth. Characteristics in PSD and SASR between four seasons for the same landforms respectively were discussed as follows.  PSD, in the sample area of the northern hemisphere, within each roughness class, all decreased in the following order: spring, summer, winter, autumn.This order was consistent with the order of Martian imbalanced ratio of four seasons (the ratio of the duration of each season). However, in the sample area of the southern hemisphere, PSD did not decrease in the same order even decreased in the disorder order (see Table 3). These characteristics were different from Earth. Theoretically. as discussed in result 2, the PSD for a day was determined by two factors: the PSD from sunrise to sunset determined by solar azimuth and the decreased PSD determined by solar elevation angle. Take the northern hemisphere on Earth as an example, since the sample area is in the northern hemisphere, the sun elevation angle is at its maximum in summer and at its minimum in winter, respectively. The sun elevation angle in spring and autumn lies between them. Besides, the sun elevation angle affects the intensity of the shielding effect (see Figure 8) The higher sun elevation angle may lead to the weaker shading effect, while the weaker shading effect causes the increase of the PSD. Meanwhile, as the sun is directly in the northern hemisphere in spring and summer, the solar azimuth may cause the total period of daytime (from sunset to sunrise) to increase, which may cause more light to be received in summer and spring. That is. the theoretical PSD without considering the shielding effect may bring certain influences to different seasons. Therefore, based on the abovementioned two reasons, the PSD in the northern hemisphere of the Earth generally decreased in the following order: summer, spring autumn and winter, which was confirmed by [12]. Similarly, we can infer that PSD in the southern hemisphere of the Earth generally decreased in the following order: winter autumn, spring and summer. However, according to the above, the seasonal characteristics of sample areas on Mars discovered in S S P was different from Earth. Then, we aimed to explore the reasons before this phenomenon. We found the reason before this finding was mainly because of the imbalanced seasons caused by Martian moving orbits and velocity (see Table 1). This unique Martian movement causes the durations of spring and summer on Mars are much longer than autumn and winter (the ratio of durations for spring, summer, autumn and winter is: 1.34: 1.25:1: 1.08, while the Earth's is nearly 1:1:1:1). Then, this imbalanced duration of seasons resulted in the weakening influences of the total shielding effect to PSDs in each season. The descending order of PSDs in the sample area of the northern hemisphere, within each roughness classes, all decreased in the consistent order with the descending order of the duration of four Martian seasons (spring, summer, winter, autumn). In the southern hemisphere, the descending order of PSDs in the sample area at first affected by the unbalanced seasons still maintains the same descending order (spring, summer, winter, autumn). Then, affected by the maximum PSD in a day for each season and the shielding effect, the descending order of PSDs gradually varied. The conclusion can be verified by the following steps: Step 1. We considered the duration of mean daily maximum PSD in each season. Two test areas were added to verify the characteristics in the southern hemisphere (the northern part and southern part of Hellas Plain). We calculated the theoretical maximum PSD in each season for these six regions and denoted as P max (added up the period from sunrise to sunset of every day in each season). Then, we divided P max in each season by its corresponding number of days to get the mean daily maximum PSDs (it is hereafter called P dmmax for brevity). In the southern hemisphere, the P dmmax in each season generally decreased in the following order: autumn, winter, spring, and summer; besides, as latitude increased, the influences of P dmmax gradually enhanced and became the leading factor (the ratio of P dmmax in four seasons ranged from 1.0002:1:1.0179:1.0177 to 1.0094:1:1.8944:1.8846). Regions in the southern hemisphere do not violate the rules mentioned above; when the sun is directly in the southern hemisphere, autumn, and winter have more daytime hours than spring and summer. It is the result of the influence of the solar azimuth.
Step 2. We considered the ratio of daily shielding effect (it may be called R S for convenience) in each season. We calculated the actual daily PSD as P a . The reduced PSD denoted as P r in each season was the difference between the Pmax and P a . Furthermore, the ratio of daily shielding effect in each season was viewed as the ratio of P r to P max . Then, we calculated the ratio of the daily shielding effect in four seasons (see Table 3). The descending order of daily shielding effect ratio was generally in the following order: summer, spring, autumn, winter. Additionally, the R S is generally less than 6%. Then, we may draw an inference: although the shielding effect has little influence on the descending order of daily PSD, the daily shielding effect has a consistent law which may lead to the same trend of seasonal PSD. It is the result of the influence of the sun elevation angle.
Step 3. Based on the above two experiments, we draw a more general conclusion: the daily PSD was influenced by theoretical daily maximum PSD and the daily shielding effect. Compared to the daily shielding effect, the theoretical daily maximum PSD was the leading factor. Along with the order of autumn, winter, spring, and summer, the mean theoretical daily maximum PSD decreases, the mean daily shielding effect leads to little influence. Then, influenced by these two determine factors, the PSD may decrease in the following trend: autumn, winter, spring, and summer. However, the imbalanced seasons impede the trend of PSD (This can be observed in Table 3, the original order of the seasons of landforms at low latitudes were the same as the order of unbalanced seasons. Then, with latitude increasing, the descending order of four seasons for corresponding landform gradually approached to the previously mentioned trend). It is the combination of these three factors that make the law of seasons different from the Earth.
Step 4. This stage was to prove the conclusions proposed in stage 3. We considered the performance of mean daily PSD in each season. We calculated the PSD of the sample regions in each season. Then, we divided the seasonal PSD of each season by its corresponding number of days to get the mean daily PSDs in each season. We found that the seasonal mean daily PSD was affected by the mean theoretical daily maximum PSD and the mean daily shielding effect. From the equator to the southern hemisphere, it quickly approached to the descending trend mentioned in stage 3: autumn, winter, spring and summer. However, due to the imbalanced seasons, the seasonal PSD at the beginning did not decrease in the same order as the seasonal mean daily PSD. It decreased in the following order: spring, summer, winter and autumn, which was the same with the descending order of imbalanced seasons. Then, with latitude increasing, it gradually varied into the sequence: autumn, winter, summer and spring, which was influenced by the mean daily maximum PSD and the shielding effect. That is, obviously, it can be concluded that although the descending order of mean daily PSDs in each season was nearly consistent with the descending order of mean daily maximum PSDs, the imbalance of seasons weakens the influences on the total PSD in each season.
Step 5. We also calculated the mean daily maximum PSDs, shielding effect ratio, and mean daily PSD for all sample areas in the northern hemisphere in each season. With latitude increasing, the influences of P dmmax gradually enhanced and became the leading factor (the ratio of P dmmax in four seasons ranged from 1.0174: 1.0176: 1: 1.0002 to 1.3419: 1.3453:1: 1.0036). The descending order of shielding effect ratio was autumn, winter, summer, and spring. Besides, the shielding effect ratio was generally less than 7%, which may lead to little influence on the descending order of total PSD in four seasons. Then, we found that under the influences of shielding effect and theoretical daily maximum PSDs, the descending order of mean daily PSDs was generally the same as theoretical daily maximum PSDs (that is, summer, spring, winter, autumn). However, due to the imbalanced Martian seasons, the descending order of PSDs in four seasons generally was: spring, summer, winter, autumn, which was in common with the descending order of imbalanced Martian seasons. That is, the inference we derived in the southern hemisphere was also confirmed in the northern hemisphere.
Note that we found that the descending order of the PSDs in S s P in four Martian seasons within each roughness class showed the same trend with the actual descending order of seasonal mean PSD no matter in the northern hemisphere or southern hemisphere. That is, S s P was able to reveal the relationship of PSD of four seasons in sample areas.
Besides, stage 2 and stage 5 provide deep support for the explanation for result 4. In the different sample areas of the northern hemisphere, along with the increase of latitude, the increase of PSD in spring and summer was generally greater than the decrease of PSD in autumn and winter. These findings lead to the result that the annual PSD decreased with latitude decreasing in the northern hemisphere.
Differently, the S s R found by us in four seasons in the sample areas showed different characteristics with the S s P. In the sample regions of the northern hemisphere, the S s R in each roughness class generally decreased in the following order: summer, spring, autumn, and winter; in the sample region of the southern hemisphere, the S s R decreased in the following order: autumn, winter, summer, and spring. Then, we considered the influences of the shielding effect and mean daily SASR in four seasons. In common with PSD, we calculated the theoretical maximum SASR as S max (the SASR from sunrise to sunset, which did not consider the shielding effect) and the daily shielding effect ratio of SASR which was similar to the method in PSD. We found that the SASR among four seasons in the northern hemisphere displayed the same descending order as that of the mean daily PSD among for season.
The reason before the phenomena that the descending order of SASR in four seasons was consistent with that of mean daily SASR in four seasons may be explained as follows. The SASR was also affected by the imbalance of seasons. But this effect may be counteracted by the shielding effect and mean daily SASR in four seasons. On the one hand, the descending order of shielding effect ratio of SASR was consistent with that of PSD. On the other hand, in each sample area, the shielding effect ratio in each season of SASR was greater than the shielding effect ratio in each season of PSD (the former minus the latter is about 0.03-1.1%). Besides, huge differences were exit on the ratio of mean daily SASR in each season, which was different from PSD (in charitum montes, the ratio of SASR in four seasons can even reach 1:1.099:5.5739:5.1697, while the ratio of PSDs in four seasons were 1.0093:1:1.8943:1.8846). That is, different from PSD, the shielding effect and mean daily SASR in each season may obviously influence the total SASR in each season. Then, under the dual effects of these two factors, the effect of the imbalanced seasons on SASR may be nearly counteracted. This conclusion was confirmed by the fact that the SASR displayed an opposite descending order of the mean daily SASR.
We extracted the decay rate of S s R and that of S s P in four seasons (the similar experimental method in result 4) and found the features similar to annual spectrums: along with latitude decreasing, the decay rate gradually decreases (see Figures 14 and 15). It may be summarized that the attenuation rate gradually decreased from the sample areas of high latitudes in the northern hemisphere to the equator, and then gradually increased from the equator to the sample areas of high latitudes in the southern hemisphere. This feature may be explained by the similar ways we provided for the same features in the annual spectrum.
Note that we may obtain a more general conclusion: SASR is more sensitive to terrain influences than PSD. This was proved by the finding that the decay rate of SASR is generally greater than that of PSD (see Figures 16 and 17). Besides, the phenomenon that the shielding effect ratio in each season of SASR was greater than the shielding effect ratio in each season of PSD in the sample area also proved this conclusion. That is, the SASR may be highly influenced by the shielding effect than PSD.
In Table 3, A represents spring, B represents summer, C represents autumn, and D represents winter. As seen in this table, under the influence of solar azimuth, the mean maximum daily PSD decreased in the following order: C, D, A, B; under the influence of solar elevation angle, the shielding effect of PSD generally decreased in the following order: B, A, C, D. Then, these two factors may become the determine factors of mean actual daily PSD. Thus, the mean actual daily PSD quickly varied to the same descending order of mean maximum daily PSD. That is, the daily PSD on Mars showed enough regularity due to the influences of solar azimuth and solar elevation angle. However, due to the imbalanced period of seasons on Mars, the seasonal PSD does not go down in the same order as the mean actual daily PSD. The seasonal PSD at the beginning decreased in the same descending order of duration in each season: A, B, D, C (the period ratio of spring, summer, winter and autumn on Mars is 1.34: 1.25:1: 1.08). Then, with the increase of latitude, the influence of sun azimuth and sun elevation angle increased; then, it caused the descending order of PSD to gradually varied. It is a gradual process. Finally, it became a fixed order: D, C, A, B. each season of PSD in the sample area also proved this conclusion. That is, the SASR may be highly influenced by the shielding effect than PSD. Figure 14. The decay rate of roughness-mean shielded astronomical solar radiation spectrums in four seasons in 12 sample areas. The decay rate of roughness-mean shielded astronomical solar radiation spectrums in four seasons in 12 sample areas. Figure 14. The decay rate of roughness-mean shielded astronomical solar radiation spectrums in four seasons in 12 sample areas. Figure 15. The decay rate of roughness-mean possible sunshine duration spectrums in four seasons in 12 sample areas. Table 3. Statistics of 6 sample areas and 2 test areas of the southern hemisphere.
The number of sample area The name of season Furthermore, the annual PSD was yielded from summing the PSDs of four seasons. That is, the above characteristics may further result in a relationship between 2s in year and 2s in the four seasons. As seen in Figure 18, the 2s in a year were generally between the maximum and minimum of 2s in four seasons in each sample area. This finding may be similar explained by how a similar manner in which we explained the descending order of 2s in the four seasons. It may be also a homogenized result of PSD in a year which generally approached to the maximum PSD in a year. However, due to the PSD in a year being a cumulative result of four PSDs of four seasons, the degree of homogenization of PSD in a year may between the maximum and minimum homogenization of PSD in four seasons.
Meanwhile, we also found the 1s in four seasons in each sample area. The increasing order of 1s in four seasons was also the same as the increasing order of shielding effect ratio. SASR was affected by the shielding effect and PSD. Then, the SASR may tend to the maximum SASR (i.e., the SASR from sunrise to sunset, regardless of shielding effect). That is, SASR may be homogenized. In common with PSD, we can infer that this homogenized characteristic may lead 1s to increase in the same order of shielding effect. Then, we proved this explanation for this finding. We calculated the variation coeffi-

Relationship between Critical Areas for Spectrums in Four Seasons
Additionally, we discovered that increasing order of X 2s in four seasons was the same as the increasing order of shielding effect ratio (see Figure 16). The reasons before this characteristic may be explained as follow. The PSD in each season generally approached to the maximum PSD in four seasons (i.e., the period from sunrise to sunset in each season, which did not consider the shielding effect). On the other hand, there may be an ideal stable space structure for PSD when PSD was not affected by the shielding effect. The weaker the shielding effect, the closer the PSD is to the maximum PSD, the higher the similarity of PSD in the grid will be, and vice versa. Afterwards, under the above, if the shielding effect enhanced, the similarity of PSDs in the grids may weaken. A larger critical area is needed to guarantee the stable spatial structure of PSD. Accordingly, the closer the PSD is to the maximum PSD, the smaller X 2s maybe. That is, the X 2s may decrease in the same order as the descending order of shielding effect.  Furthermore, the annual PSD was yielded from summing the PSDs of four seasons. That is, the above characteristics may further result in a relationship between X 2s in year and X 2s in the four seasons. As seen in Figure 18, the X 2s in a year were generally between the maximum and minimum of X 2s in four seasons in each sample area. This finding may be similar explained by how a similar manner in which we explained the descending order of X 2s in the four seasons. It may be also a homogenized result of PSD in a year which generally approached to the maximum PSD in a year. However, due to the PSD in a year being a cumulative result of four PSDs of four seasons, the degree of homogenization of PSD in a year may between the maximum and minimum homogenization of PSD in four seasons.
classes: 0-5 m, 5-10 m, 5-10 m, 10-15 m, 15-20 m, 20-25 m ,25-30 m, 35-40 m, 40-45 m, 45-50 m ,50-55 m, 55-60 m, 60-65 m, 65-70 m, 70-75 m, 75-80 m, 80 m-) In E 21, an 10 m equal-interval roughness classification scheme is adopted (i.e., the roughness is classified into 28 classes: 0-5 m, 5-10 m, 5-10 m, 10-15 m, 15-20 m, 20-25 m ,25-30 m, 35-40 m, 40-45 m, 45-50 m ,50-55 m, 55-60 m, 60-65 m, 65-70 m, 70-75 m, 75-80 m, ,80-85 m, 85-90 m, 90-95 m, 95-100 m, 100 m-). As presented in the figure, when the classified roughness is greater than 80 m (lie in the 18,19,20,21 class), the spectrum showed a tendency to be unstable. This is mainly because when the critical region is expanding, few grids within the roughness class of 18-21 are obtained (as observed, only about 1-3 number of the corresponding grids in each class were obtained). These grids can be regarded as singular values and are not representativeness of statistical significance.  4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27  Meanwhile, we also found the X 1s in four seasons in each sample area. The increasing order of X 1s in four seasons was also the same as the increasing order of shielding effect ratio. SASR was affected by the shielding effect and PSD. Then, the SASR may tend to the maximum SASR (i.e., the SASR from sunrise to sunset, regardless of shielding effect). That is, SASR may be homogenized. In common with PSD, we can infer that this homogenized characteristic may lead X 1s to increase in the same order of shielding effect. Then, we proved this explanation for this finding. We calculated the variation coefficient of SASRs in each season in each sample area. We found that the descending order of the variation coefficient of SASR in each sample area is the descending order of the shielding coefficient. That is, the larger the shielding coefficient, the more complex the distribution of SASR in the sample area. The complex distribution of SASR may lead to a larger X 1s .
Regardless of the year or season, the value of X 1s was generally greater than X 2s (the ratio of X 1s to X 2s ranged from 1.12 to 2.172). This finding can be explained as follows. Based on Equation (13), the PSD is a principal parameter of SASR which may greatly influence SASR. Besides, since there are other involved parameters in the calculation of SASR, (such as slope, aspect.), the distribution of SASR may be more complex than PSD. When we found the X 2s , due to the more complex SASR, a larger area may be required to obtain the similarities in the distribution of SASR. Accordingly, for the same duration (a year or each season), X 1s may be generally greater than X 2s .
Furthermore, in common with X 2s , the X 1s in a year was between the maximum and minimum of X 1s in four seasons. The reason behind this characteristic of X 1s may be explained in the same way we explained the reason for the similar characteristic in X 2s . These findings were confirmed by the following Equation: X 2s−year = X 2s−spring +X 2s−summer +X 2s−antumn X 2s−winter )/4 R 2 = 0.74 (41) X 1s−year = X 1s−spring +X 1s−summer +X 1s−antumn X 1s−winter /4 R 2 = 0.79 (42)

Discussion
Combining planetology, Martian solar radiation theory, and the astronomical radiation theory on Earth, a complete process to calculate the Martian SASR and PSD was proposed in this paper. Compared to the previous research on Mars, this method took the terrain influences into account for the first time. That is, to some extent, it improved the accuracy of the previous models. Thus, this method is promising for future works about SASR on Mars. The spectrum using the slope as the x-axis was not suitable to be adopted on Mars. Thus, this paper proposed the new concept of spectrums for PSD and SASR, which was tested to be suitable for the study of SASR and PSD on Mars. The parameter to describe terrain relief is roughness. The RMS, which was generally viewed as a terrain parameter to describe roughness, was adopted to reveal the roughness on the Martian surface. However, although RMS was adopted as a terrain parameter in common with the slope to describe planetary surface topographic relief, it has never been used in spectral studies. According to the finding of the result, it is clear that S s R and S S P can be efficient in exploring the special spatial-temporal distribution of SASR and PSD. The discussion about findings was as following respects: (1) SASR and PSD showed the law of gradual variation with terrain relief. Either on an annual scale or a quarterly scale, the value of S s R and S s P in different Martian landforms generally decreased. This phenomenon may reveal that SASR and PSD were highly influenced by terrain relief. In each sample area, the rugged terrain may cause the decrease in SASR and PSD. S s R and S S P in a year or four seasons revealed the interaction of SASR and PSD with terrain relief.
(2) Under the influences of shielding effect, the seasonal and annual spatial structure of SASR or PSD on Mars revealed by S s R and S s P showed the latitude anisotropy characteristics. The descending order of S s R and S s P in four seasons showed a certain regularity: the attenuation rate decreased from the sample areas of high latitudes in the northern hemisphere to the equator, and then increased from the equator to the sample areas of high latitudes in the southern hemisphere. This regularity is affected by the sun elevation angle and sun azimuth, which are further affected by the shielding effect. That is, it is still a manifestation of the terrain influences on SASR and PSD.
(4) As discussed in the results, according to the findings of S s R and S s P, the law of the SASR and the PSD on Mars varied with latitude are totally different from those on Earth, either on an annual scale or on a quarterly scale. This was mainly because of the imbalanced seasons of Mars. Besides, SASR is more regular than PSD. The reason before this feature may be because SASR is more sensitive to the shielding effect than PSD. In all sample areas, the finding that attenuation rate of S s R was greater than S s P and the finding that the shielding effect ratio of SASR was greater than PSD provided deep support for this explanation.
(5) The X 1s and X 2s are effective indexes to reflect the stable spatial structure of SASR and PSD respectively in different geomorphology. The X 1s and X 2s may be regarded as minimum test areas in this sample area for the research on SASR or PSD since it showed sufficient stability when the test area is larger than those. Besides, the X 1s and X 2s in a year seemed to be positively correlated with the variable coefficients of roughness. And the X 1s and X 2s in a year may be the mean of that in four seasons. The increasing order of X 1s and X 2s in four seasons was the same as the increasing order of shielding effect ratio. These relations may help us quickly determine the minimum test area of the sample area.
Note that the terrain analysis on Mars provides additional insights beyond planetary, astronomical and climatological parameters (which are not within the scope of this issue).

Discussion on the Influence of Roughness Classification Scheme on the Spectrum
The previous studies about spectrum slope classification scheme generally based on the following two points: in the spectrums proposed by previous studies, equal-interval slope classification was generally used as a classification scheme better than other schemes to ensure that comparative studies can be conducted in each landform; besides, 45-90 • were generally viewed as a single class since the grids with the value of slope between 45-90 • accounts for little of the total grid area [11,12,78,79,133].
Referring to the above spectral design, we calculated the roughness of 50 randomly selected test areas and 12 sample areas on Mars and found the number of grids with the value of RMS greater than 80 m accounted for less than 6% of all grids (ranged from 0% to <5.83%) in any test area and sample area. Thus, for the roughness spectrum, we design the structure of roughness classification as a 5 m equal-interval roughness classification scheme the spectrums whose x-axis was slope, since the number of grids whose RMS value was greater than 80 m accounted for less than 6% of all grids in any sample area, RMS whose value was greater than 80 m was used as a single class in the spectral scheme). We also seek to further classify the grids whose RMS is greater than 80 m. It was not able to be extracted in some regions of flat terrain such as impact craters or planitia (because there are no grids whose RMS was greater than 80 m) or display the disorder trend for different S s R and S s P in the class whose RMS value was greater than 80 m (see Figure 17).
We discovered the corresponding S s R and S s P in each sample area by adopting the 5 m equal-interval roughness classification scheme (called E 17 for convenience, see Figure 17). Then we seek to explore different influences of different roughness classification on the spectrums. we also investigated the corresponding S 1s P and S 2s P in each sample area by adopting the 3 m equal-interval roughness classification scheme (called R 28 for convenience, see Figure 18) and 10 m equal-interval roughness classification scheme (called E 9 for convenience, see Figure 18).
The corresponding characteristics in E 17 can be also found in the E 28 and E 9: e.g., along with the roughness class increasing, the corresponding S s R and S s P in a year and in four seasons generally decreased. The relationship of S s R and S s P in sample area in four seasons based on E 24 and E 9 were also consistent with the law of E 17. The explanation for this result can be in the same way as E 17.
The same spectral characteristics found on the corresponding spectrums based on different equal-interval roughness classification reveals the reliability of the characteristics in SASR and PSD which we found on Mars. In other words, it seems that the spectrum methods proposed by us show a certain remarkable ability in accurately revealing the spatial-temporal distribution characteristics of SASR and PSD.

The Commonalities and Differences between the Two Spectrums and Other Spectrums Proposed Before
We discussed the commonalities and differences of the proposed spectrum method and other spectrums before as follows.
The spectrum method generally adopted a form of two-dimensional spectrum to reveal the stable spatial structure of specific matter composition (slope, PSD or SASR). Namely, they usually performed as a group of the histogram.
The morphology of previous spectrum method, S S R and S S P is closely related to the extracted region. For each style of spectrum, the critical area generally exists in different landforms. Each spectrum can only be found when the extracted region was larger than the critical area. Besides, when the extracted region was larger than the critical area, the spectrums extracted in any location of sample area generally present self-similarity and stability. Furthermore, the slope spectrum is the theory of geomorphological principles, which is employed to discover the stable structure of slope. The innovation of the other derivative spectrums is the application of geomorphological principles to astronomy, which were proposed to explore the stable spatial structure of SASR or PSD.
(2) In Findings The spectrums may be viewed as the representative features of different landforms either on Mars or Earth. The spectrum method was successfully employed to explore the stable spatial structure of SASR or PSD in different landforms. The revealed stable structure of different landforms was generally closely related to its topography.
The spectrums revealed the temporal distributions of SASR or PSD for the corresponding landform. For each landform, as observed in the spectrums, the spatial structure of SASR or PSD in different temporal scale is different.
In the spectrums either on Mars or Earth, the SASR or PSD generally showed a downward trend in the corresponding spatial structure which revealed the interaction of SASR or PSD with terrain influences.
Besides, the SASR and PSD revealed by the spectrums generally showed a highly correlations with the shielding effect caused by the rugged terrain in different landforms.
The critical area of the spectrums seemed to be positively correlated to terrain complexity. Furthermore, certain relations seemed to be found between the critical areas of different temporal scale.

Differences (1) In Characteristics
The application scope is different. The other spectrums cannot be found on Mars. The S s R and S s P was firstly proposed to be employed in the medium resolution DEM data for conducting the terrain analysis, which is convenient for the SASR and PSD study in Mars. The other spectrum was used to be employed in the high-resolution DEM data in Earth, which is not suitable for Mars.
The employed methods and application aim of the spectrums is different. The slope spectrum is the theory of geomorphological principles, which is employed to discover the stable spatial structure of the slope. The other spectrum originated from slope spectrum to explore the SASR or PSD. The slope-mean SASR spectrum and slope-mean PSD ratio spectrum adopted slope to explore the stable spatial structure of SASR or PSD shielding effect under terrain relief. According to the current situation of Martian study, the S s R and S s P were innovatively proposed to use roughness to describe the spatial structure of SASR and PSD under terrain relief. The RMS, which was generally regarded as a terrain parameter to describe roughness, was adopted to reveal the roughness on Martian surface. However, although RMS was usually adopted as a terrain parameter in common with slope to describe topographic relief, it has never been used in spectral studies.
Thus, the composition of the spectrum is different. Slope spectrum performed as a group of the histogram, whose x-axis is slope, y-axis is slope percentage. Slope-mean SASR spectrum performed as a group of the histogram, whose x-axis is slope, y-axis is mean SASR; slope-mean PSD ratio spectrum is a group of the histogram whose x-axis is slope, y-axis is the ratio of the mean PSD considering terrain influences on the mean PSD on the horizontal plane; S s R and S s P performed as a group of the histogram, whose x-axis is RMS, y-axis is mean SASR or PSD.
(2) In Findings The spatial-temporal variations and laws of Mar revealed by spectrums was generally more complicated than Earth.
Different from the features of SASR and PSD on Earth revealed by [11,12], the SASR and PSD on Mars showed the unregular latitude anisotropy characteristics found in part 3. The descending order of SASR and PSD in four seasons was also different form Earth. As discussed in the finding of part 3, the regularity is mainly because the imbalanced seasons caused by the differences of velocity and rotations between Mars and Earth.
In Mars, the spatial structure of SASR seemed to be more sensitive to the shielding effect than PSD unlike Earth, which is proven by our result.
The critical area of SASR and PSD was generally greater than that of [11,12].The main reason for this is the total area of Martian landforms is generally larger than that of Earth.
The previous spectrums which based on slope area generally unable to be extracted on Mars. The spectrums based on RMS is suitable for Mars.

Limitations
The landform which has too flat terrain may lead to the failure of spectral extraction. We failed to extract the spectrum of dune fields though we have tried all approval dune fields (Abalos Undae, Aspledon Undae, Hyperboreae Undae, Ogygis Undae, Olympia Undae, Siton Undae). The mean value of the roughness for grids was less than 5 m. The percentages of the part of the area with roughness within the interval of 75-80 were between the value of 0 and 0.0015% and that of larger than 80 m were between the value of 0 and 0.00145%. Besides, the percentages of the part of the area with roughness larger than 30 m were less than 0.003% in the six dune fields. Obviously, the dune fields generally have a flat terrain. Then, we found that the roughness classes of the area which account for low percentages were failed to derive the corresponding value. Namely, when the topographic feature information of the landform is too little, we cannot derive the corresponding spectrums. On the other hand, it reflects the close relationship between the spectrum and the geomorphology.

Application of this Study
The models of SASR on Mars proposed in this paper can be used for future studies of astronomical radiation on Mars. Previous methods focus on the impact of the Martian meteorological environment on astronomical radiation but ignored the consideration of the terrain influence. This paper took the terrain influence into consideration, which is applicable for future research on Mars.
X 1s and X 2s may help researchers to determine the minimum test region of the sample areas. The spatial structure of SASR and PSD are stable only when the extracted region of the sample area is larger than the corresponding spectral critical area. Thus, they may be able to be as minimum test areas for corresponding research on SASR and PSD. If the extracted regions used to study SASR and PSD are less than the corresponding X 1s and X 2s , the obtained spatial structure of SASR and PSD may not reflect the distribution characteristics of those in this area. Moreover, in the 12 sample areas, sufficient positive correlation generally exists between X 1s and X 2s in a year and terrain relief (see Section 3.1). Thus, when we conducted research on the spatial distribution of SASR or PSD in a year, using the Equations proposed in 3.1 and the variable coefficient of roughness, we may be able to obtain minimum test regions for sample areas without the complex calculation process of obtaining the X 1s and X 2s S s R and S s P, developed from the spectrum method in geoscience, were appropriate methods to reveal the spatial-temporal distribution of SASR and PSD on Mars. Different from the statistic method and the comparative method, S s R and S s P can reveal the complicated interaction between terrain relief and SASR or PSD. S s R and S s P in different areas generally showed different features, i.e., they respectively were able to reveal the differences of SASR or PSD between different regions under terrain relief.
The value of SASR in S s R in each RMS class is the characteristic of the stable spatial structure of the corresponding sample area. That is, it can be regarded as the SASR's characteristic value in each RMS class of the corresponding sample area. Similarly, the value of PSD in S s P in each RMS class can be regarded as the PSD's characteristic value in each RMS class of the corresponding sample area. Through the combination of roughness and S s R and S s P, we can quickly obtain the corresponding value of SASR and PSD in the sample area. The value of astronomical radiation and possible sunshine duration in a year or four seasons may be quickly calculated by S s R and S s P proposed by us. The calculations for astronomical radiation and possible sunshine duration were so complex that it often requires the application of computer parallels to improve efficiency.

Conclusions
The models of SASR and PSD that can be applied to Mars are proposed. They compensate for the shortcomings of SASR and PSD models in previous Martian studies, taking into account terrain influences and improving time accuracy. These models can further extend to other Martian solar radiation models.
Different from the previous traditional geographic method such as the statistic method and comparative method, an innovative spectral method was proposed and applied to the study of SASR and PSD. We conducted experiments to extract the S s R and S s P of 12 Martian landforms. Then, we draw the following conclusions: (1) The spectral method is a quantitative method to be more effective in identifying and characterizing the spatial-temporal distribution of SASR or PSD in sample areas. The seasonal combination spectrum in four seasons is an effective qualitative description of the temporal distribution of SASR and PSD in different landforms. (2) S s R and S s P revealed the complex interactions between the SASR or PSD and the terrain relief. Under the terrain influences, in the S s R and S s P of the same landforms, the SASR and PSD showed a downward trend. This feature revealed that the SASR and PSD tended to decay under the influences of shielding effect caused by terrain relief. (3) SASR and PSD showed the latitude anisotropy characteristics. The latitude anisotropy characteristics discovered on Mars were complex and different from Earth due to the imbalanced seasons. In essence, this feature is also a manifestation of different shading effects caused by solar elevation angle. (4) SASR is more sensitive to the shielding effect than PSD which is proved by the corresponding experiments. Based on it, SASR showed more regular laws than PSD under terrain relief in a year or four seasons. (5) X 1s or X 2s can be a parameter to determine the minimum test regions for SASR or PSD of sample areas. The spatial structure of SASR or PSD become stable if the extracted region was larger than X 1s or X 2s . The relations discovered in results may help us to quickly found and test them.
Since the slope may be not appropriate to reveal the terrain relief on Mars, the method to extract X 1s and X 2s adopted roughness to successfully reveal the relations between terrain relief and SASR or PSD. This new spectral method theory sheds new light on the spatial-temporal structure of SASR and PSD and may provide a reference for the further future work of SASR and PSD on Mars.
The corresponding experiment result herein may prove that this method may be beneficial for us to explore the spatial-temporal distribution of SASR and PSD in different landforms on Mars. This method lays the groundwork and can be further extended to the study of the spatial-temporal structure of other solar radiation types.
The application of spectral method in this paper expands the application fields and techniques of spectral method. Additionally, our work advances the modeling and theoretical analysis of the solar radiation of Mars. It makes up for the deficiencies in the current relevant models and analyses.
Author Contributions: Conceptualization, Siwei Lin and Nan Chen; methodology, Siwei Lin and Nan Chen; software, Siwei Lin; validation, Siwei Lin and Nan Chen; formal analysis, Siwei Lin; investigation, Siwei Lin and Nan Chen; resources, Siwei Lin and Nan Chen; data curation, Siwei Lin and Nan Chen; writing-original draft preparation, Siwei Lin and Nan Chen; writing-review and editing, Siwei Lin and Nan Chen; visualization, Siwei Lin and Nan Chen; supervision, Nan Chen; project administration, Nan Chen; funding acquisition, Nan Chen. All authors have read and agreed to the published version of the manuscript.
Funding: This work was supported by the National Natural Science Foundation of China (grant numbers 41771423, 41491339, 41930102, and 41601408) and by the industry-university-research cooperation project for the social development of Fujian Province, China (grant number 2018Y0054).