Constructing a High-Accuracy Geometric Model for Moon-Based Earth Observation

: The Moon provides us with a long-term, stable, and unique place for Earth observation. Space agencies of various countries, including the United States, China, and Italy, have made the realization of Moon-based Earth observation an objective of their lunar missions. To date, although some conceptual studies have been presented, an accurate geometric model for Moon-based Earth observation has not yet been described. This paper introduces such a geometric model, which connects the attitude of a Moon-based sensor with a corresponding ﬁeld of view on the Earth’s surface. The aberration and light time correction are involved. Due to the lack of high-quality experimental data, one qualitative comparison with Chang’E lander images and another quantitative comparison with software output are made. The comparison results show good similarity. The overall model accuracy is evaluated to be better than 4 (cid:48)(cid:48) at current stage and will be better than 1.5 (cid:48)(cid:48) if the seleno-graphic position is accurately determined. In direct geolocation process, the aberration and light time will cause a total 0.7 km deviation on the ground near the sublunar point. In SAR range history simulation, the light time e ﬀ ect will lead to a linear error, as large as tens of meters, throughout the integration time.


Introduction
No matter how far humans want to travel in the universe, the Moon, the only natural satellite of the Earth and our nearest neighbor, is the first place we need to colonize. Since 1958, satellites and spacecrafts have visited the Moon more than 100 times. Countries and institutions including the Soviet Union-Russia, USA, ESA, China, Japan, and India have been involved in these explorations. In this century, more than 10 missions have been successfully launched, and tens of missions are scheduled. With technological progress and increases in our understanding of the Moon, the colonization of the Moon is no longer an envisaged possibility but is a real challenge humans face. Recently, the ESA proposed a "moon village" program, which is a Moon base under international collaboration among spacefaring nations for science, business, mining, and even tourism [1]. NASA also started its new lunar mission, named Artemis, to return astronauts to the lunar surface by 2024 [2]. Similar plans have emerged over past few decades, but until now, humans have not truly had the ability to realize these plans. Such a permanent settlement on the lunar surface will be significant to scientific studies and deep space exploration. For Earth science, it will also provide an opportunity to complete our Earth In addition, the meanings of several other coordinate systems have to be specified. Similar to the geographic coordinate system (GCS), the seleno-graphic coordinate system (SGS) is set up to define Lat-Lon positions on the lunar surface, and its origin is considered the intersection point of the X-axis of the body-fixed coordinate system and the lunar reference sphere. The lunar reference ellipsoid is a sphere with a radius of 1737.4 km. In addition, the three other right-handed three-dimensional space coordinate systems used in the following sections are: Sensor Coordinate System (SCS): The SCS is a sensor-fixed reference frame. The origin is the optical center for the camera or the phase center for the antenna. The Z-axis is the optical axis or the direction of the main lobe, namely, the line of sight (LOS) in the following sections. The XY plane is perpendicular to Z. The directions of the X-axis and Y-axis refer to the antenna shape or the image coordinates. In this paper, they are implicitly defined by the zenith angle and azimuth angle.
Observatory Local Coordinate System (OLCS): The origin of the OLCS is an artificially defined center of the observatory with an elevation of H. The Z-axis is the same as the local normal vector, while the XY plane is the local tangent plane. The X-axis points to the north. The zenith angle-and the azimuth angle-related SCS and OLCS will be introduced in Section 3.4.

Selenocentric Body-fixed Coordinate System (SBCS):
The coordinate origin is the lunar mass center. Principle axis (PA) coordinates and mean-Earth/rotation-axis (MER) coordinates, which are related by a constant three-angle rotation, are two commonly used body-fixed coordinate systems [23]. The PA coordinates are directly linked with the ICRS by physical librations [24], while lunar cartography customarily refers to MER [25].

Considerations Concerning the Time Coordinates
The Earth orientation parameter (EOP) data and JPL ephemeris used in this paper are presented in different time coordinates. The terrestrial time (TT) coordinating the Earth's attitude has a mean rate close to the proper time of the rotating Earth and can accurately be represented by International Atomic Time (TAI), an idealized atomic time system with the SI second on the geoid as the unit of measurement [26]: TT = TAI + 32.184 second (1) Barycentric dynamical time (TDB) serves as an independent time argument of barycentric ephemerides, such as the JPL ephemeris. The relationship between TDB and TT is given by [27]: where P is a periodic term with a main component with an approximately 1.7 ms amplitude and a one-year period. P leads to periodic geometry deviation during the Earth's orbital cycles; however, its amplitude is too small for the Earth observation application. The maximum space deviation of the lunar mass center is less than 2 m, much smaller than other errors analyzed by Section 5.3. Universal Time Coordinated (UTC) system, which is the standard for all civil time systems, is different from TAI because the former has integer leap seconds forecasted by the International Earth Rotation Service (IERS). The time conversion is realized by the Standards of Fundamental Astronomy (SOFA) time scale and calendar tools in this paper.

JPL Development Ephemeris and Earth's Attitude
The JPL ephemeris is in good agreement with lunar laser ranging (LLR), which now achieves millimeter-scale accuracy [28]. In this paper, we use DE430 to calculate the positions and velocities of the helio-center, geocenter, and seleno-center in the ICRS, as well as the amplitude and derivative of the libration. The position vector should be translated from the ICRS to the GCRS for the next step.
Usually, the position of the reflector on the Earth is given in the GCS, defined on the Earth's reference ellipsoid. The coordinates in the ITRS can be calculated by using the geographic coordinates and parameters of the ellipsoid. Since the center of the Sun, Earth, and Moon are given in the GCRS, Remote Sens. 2019, 11, 2611 5 of 25 we need to understand the orientation matrix between the ITRS and GCRS that describes the Earth's attitude in the universe. According to IERS conventions, the orientation matrix M ITRS_GCRS is presented by [29]: M ITRS_GCRS = M P−N ·M ER ·M PM (3) where the M P−N is the transform matrix of the precession and nutation, M ER is the transform matrix of the Earth's rotation, and M PM is the transform matrix of the polar motion. The matrix M P−N ·M ER ·M PM transforms the position vector from the ITRS to the GCRS. The IERS conventions shows how to use the models and data to calculate M ITRS_GCRS . The SOFA tools for Earth attitude provide the program that implements the official IAU algorithms based on EOP data.

Orientation Matrices
The relation between the sensor's position and pointing direction is given by a series of transformation and an intersection. The flow chart of simulation is presented in Section 4. The orientation matrices for transformation are introduced in this section, and the process of transformation will be introduced in the next section. Assuming the coordinate system rotates θ clockwise around the X, Y, and Z axes, the rotation matrices are: The vector in the new coordinate system can be calculated by pre-multiplying the unit rotation matrix to the vector in the original coordinates system.
For clarity, we assumed the sensor is on a two-dimensional altazimuth turntable [30,31]. In its initial state, the SCS and OLCS are the same. When controlling the direction of the LOS, we first rotate the azimuth angle α I clockwise around the Z-axis and then rotate the zenith angle θ I clockwise around the Y-axis. In some cases, a third rotation is needed. For example, if we use a rectangular antenna, the long side of the footprint is expected to be parallel with the meridian, and an extra T z should be added. The rotation angle −α R is the clockwise angle around the Z-axis, and it is zero for circular aperture. Thus, the orientation matrix M SCS_OLCS (Figure 1) can be written as: If the observatory position is (ξ M , ∅ M ) in the SGS, we can get the orientation matrix M OLCS_MER ( Figure 1).
The orientation matrix M MER_PA connecting the MER to PA is: The angles β 1 , β 2 , and β 3 refer to the constant parts of the three libration parameters. Their values compatible with DE430 are 67.573 , 78.580 , and 0.285 [24]. If the position of the observatory is directly given by PA coordinates x PA , y PA , z PA T , such as the coordinates of the LLR retroreflectors, the orientation matrix from the OLCS to the PA can be directly expressed by substituting the PA latitude and longitude for the MER latitude and longitude in Equation (6).
The three Euler angles of libration areΩ L ,î L , andû L . On the celestial sphere, let N 1 stand for the ascending node of the lunar mantle equator to the ICRS equator ( Figure 2). Thus, it can be inferred that the three Euler angles are the angle from the X-axis of the inertial frame along the ICRS equator to the intersection of the lunar mantle equator, the inclination of the lunar mantle equator from the ICRS Remote Sens. 2019, 11, 2611 6 of 25 equator, and the longitude from N 1 to the lunar prime meridian along the lunar mantle equator [25,32]. Thus, the orientation matrix M PA_GCRS can be presented by: that describes the Earth's attitude in the universe. According to IERS conventions, the orientation matrix _ is presented by [29]: where the is the transform matrix of the precession and nutation, is the transform matrix of the Earth's rotation, and is the transform matrix of the polar motion.    The orientation matrix _ connecting the MER to PA is: The angles , , and refer to the constant parts of the three libration parameters. The orientation matrix M GCRS_ITRS is already given in Section 3.4 and is calculated by SOFA. Note that the orientation matrix is orthogonal and has the following property:

Transformation from the SCS to the GCRS
In the SCS, the LOS can be presented by: In the OLCS, the LOS can be presented by: The position of the sensor center (PSC) in the OLCS is − −−− → S OLCS . This vector relies upon the system design. For simplicity, we assume the two centers coincide in the simulation.
In the MER, the LOS can be presented by: And the PSC can be presented by: where the R M is the radius of the Moon, and H is the elevation of the observatory center.
In the PA, the LOS and PSC can be presented by: In the GCRS, they are presented by: where − −−− → M GCRS is the position of the Moon in the GCRS.

The Light Aberration
Light takes approximately 1.3 s to travel from the Earth to the Moon, which is much longer than the spaceborne observation. Thus, an obvious geolocation problem for Moon-based Earth observation is that the geometry changes from the time when the light leaves the Earth to the time when the light arrives at the sensor. This phenomenon, known as light aberration in astronomy, results from the relative motion between the reflector and the observer. Recently, high-resolution remote sensing processing has taken aberration compensation into consideration [33,34]. This effect is calibrated by introducing a corrective rotation. Intrinsically and more precisely, a two-step method is developed for moon-based earth observation, and it also can be used in spaceborne cases [35]. For Earth observation, the aberration can be divided into two parts. The part related to the motion of the sensor will be Remote Sens. 2019, 11, 2611 8 of 25 discussed in this section, and the part caused by the Earth's rotation, called the light time effect, will be discussed in Section 3.7. The aberration angle can be approximately presented by: where θ S is the viewing angle related to the speed vector, c is the speed of light, and v is the speed of the observer. If the observing angle is perpendicular to the relative speed, the angular aberration can easily be derived and is v/c. For example, if the orbital speed of the Moon is 1.02 km/s, the viewing angle is 90 • , and the distance between the reflector and the seleno-center is 380,000 km; thus, the light aberration is approximately 3.4 × 10 −6 rad (0.7 ), corresponding to a 1.3 km deviation on the Earth's surface.
In the GCRS, the velocity of the sensor is v . Therefore, the speed addition can be described by the Lorentz transformation: where v is the magnitude of where − −−−−−− → VM GCRS is the velocity of the seleno-center in the GCRS, and: Formulating Equation (22) is a complex task [36,37]. For some applications, dΩ L dt and dî L dt can be considered approximately 0, and dû L dt can be considered 13.2 degrees per day [10]. In this paper, the result is numerically computed by using the Chebyshev polynomial to interpolate the JPL ephemeris.

The Intersection Point
In the previous calculation, the time-variant parameters refer to the receiving time t 2 , so t 2 is not explicitly written out in the equations. However, the position of the reflector refers to the 'departure time' t 1 , namely, the time when light leaves the reflector, in the geometry. Thus, the time difference should be clearly represented in the following equations. The departure time can hardly be measured by an Earth observation system, but the relationship below stands ( Figure 3): where L = c(t 2 − t 1 ) is the length of the LOS. 10 is not explicitly written out in the equations. However, the position of the reflector refers to the 'departure time' , namely, the time when light leaves the reflector, in the geometry. Thus, the time difference should be clearly represented in the following equations. The departure time can hardly be measured by an Earth observation system, but the relationship below stands ( Figure 3): where = ( − ) is the length of the LOS.

S(t 1 ) S(t 2 )
A(t 1 ) Meanwhile, at in the GCRS, ( ) ⃗ is on the terrestrial ellipsoid, so we get: simple approximation with a max error less than 10 m is presented. Note that if the LOS and the ellipsoid intersect, the range of is less than 21.5 ms because: where is the length of ( ) ⃗ . We can estimate by replacing  Meanwhile, at t 1 in the GCRS, is on the terrestrial ellipsoid, so we get: ) T of the intersection point of the LOS and the terrestrial reference ellipsoid in the GCRS at time t 1 , and L = c(t 2 − t 1 ) is the length of the LOS. The major semi-axis a e is 6378.137 km and the minor semi-axis b e is 6356.752 km in terms of WGS84 ellipsoid. By solving the equation, we can determine the position of the reflector (24) only has numerical solutions. Here, a simple approximation with a max error less than 10 m is presented. Note that if the LOS and the ellipsoid intersect, the range of t 1 is less than 21.5 ms because: where s is the length of . We can estimate t 1 by replacing M GCRS_ITRS (t 1 ) with in Equation (24). In this time window, the ITRF is under-or over-rotated with a 10 m max error (the Earth's mean rotation speed is approximately 7.29 × 10 −5 rad/s) corresponding to an error of 3.3 × 10 −8 s in the estimation of t 1 . In this paper, is chosen because the case wherein the LOS is pointed to the geocenter is mainly discussed.
In particular, imaging spectrometer models are usually described by the colinear equation, which is an effective planar approximation within limits. To address the whole disk, we follow the same procedure but substitute (10): where f is the focal length. Then, all image coordinates (x, y) are connected to geographic coordinates. This method is computationally inefficient but has a high accuracy.

Geometry Simulation
The simulation procedure follows the model specified from Sections 3.4-3.7. DE430 from JPL and EOP 14 C04 from IERS are two input data tables. Usually, times are initially presented in UTC and should be converted to TDB and TT to interpolate the two input tables. The calculated results are tested by the attached test data and the results from the IERS Rapid Service/Prediction Center (http://maia.usno.navy.mil/t2crequest/t2crequest.html). Due to the libration and the ellipticity of the orbit of the Moon, the LOS will shift out of the Earth disk if the azimuth angle and zenith angle do not change with the revolution of the Moon. To make Equation (24) always solvable, a simple strategy is to point the LOS to the geocenter: Then, the coordinates of the intersection point and corresponding azimuth angle and zenith angle can be calculated by a slight adjustment of the proposed method ( Figure 4). The coordinates of the lunar surface are directly given in the PA. The following analyses use these criteria, unless otherwise specified.
to the geocenter is mainly discussed.
In particular, imaging spectrometer models are usually described by the colinear equation, which is an effective planar approximation within limits. To address the whole disk, we follow the same procedure but substitute _ ⃗ for ⃗ in Equation (10): where is the focal length. Then, all image coordinates ( , ) are connected to geographic coordinates. This method is computationally inefficient but has a high accuracy. The simulation procedure follows the model specified from section 3.4 to 3.7. DE430 from JPL and EOP 14 C04 from IERS are two input data tables. Usually, times are initially presented in UTC and should be converted to TDB and TT to interpolate the two input tables. The A lack of experimental data makes it very hard to validate simulation results. Optically imaging the Earth is not one of the objectives of the Chang'E-3 mission. Although the Terrain Camera (TCAM) lander did take several photos of the Earth, the spatial resolution is only approximately 70 km. Since these images are the only measured data accessible to us, we compromise with a comparison of the imaged Earth disk and the simulated Earth disk ( Table 1). The elevation of the lander is −2632.0 m, and the location is (44.1206, −19.5124) in MER.DE421 [38], i.e., (44.1412, −19.5382) in PA.DE430. The orientation angles differ slightly from the angles in Equation (7) [32]. In addition, the simulation results from the STK 11 are introduced as mutual authentication data. A more credible validation method will only rely on field measurements in future missions. calculated results are tested by the attached test data and the results from the IERS Rapid Service/Prediction Center (http://maia.usno.navy.mil/t2crequest/t2crequest.html). Due to the libration and the ellipticity of the orbit of the Moon, the LOS will shift out of the Earth disk if the azimuth angle and zenith angle do not change with the revolution of the Moon. To make Equation (24) always solvable, a simple strategy is to point the LOS to the geocenter:

Geometry Simulation
Then, the coordinates of the intersection point and corresponding azimuth angle and zenith angle can be calculated by a slight adjustment of the proposed method ( Figure 4). The coordinates of the lunar surface are directly given in the PA. The following analyses use these criteria, unless otherwise specified.
Then, the coordinates of the intersection point and corresponding azimuth angle and zenith angle can be calculated by a slight adjustment of the proposed method ( Figure 4). The coordinates of the lunar surface are directly given in the PA. The following analyses use these criteria, unless otherwise specified.
Then, the coordinates of the intersection point and corresponding azimuth angle and zenith angle can be calculated by a slight adjustment of the proposed method ( Figure 4). The coordinates of the lunar surface are directly given in the PA. The following analyses use these criteria, unless otherwise specified.
Then, the coordinates of the intersection point and corresponding azimuth angle and zenith angle can be calculated by a slight adjustment of the proposed method ( Figure 4). The coordinates of the lunar surface are directly given in the PA. The following analyses use these criteria, unless otherwise specified.
Then, the coordinates of the intersection point and corresponding azimuth angle and zenith angle can be calculated by a slight adjustment of the proposed method ( Figure 4). The coordinates of the lunar surface are directly given in the PA. The following analyses use these criteria, unless otherwise specified.
Then, the coordinates of the intersection point and corresponding azimuth angle and zenith angle can be calculated by a slight adjustment of the proposed method ( Figure 4). The coordinates of the lunar surface are directly given in the PA. The following analyses use these criteria, unless otherwise specified.
Then, the coordinates of the intersection point and corresponding azimuth angle and zenith angle can be calculated by a slight adjustment of the proposed method ( Figure 4). The coordinates of the lunar surface are directly given in the PA. The following analyses use these criteria, unless otherwise specified.
Then, the coordinates of the intersection point and corresponding azimuth angle and zenith angle can be calculated by a slight adjustment of the proposed method ( Figure 4). The coordinates of the lunar surface are directly given in the PA. The following analyses use these criteria, unless otherwise specified.
Then, the coordinates of the intersection point and corresponding azimuth angle and zenith angle can be calculated by a slight adjustment of the proposed method ( Figure 4). The coordinates of the lunar surface are directly given in the PA. The following analyses use these criteria, unless otherwise specified.
Then, the coordinates of the intersection point and corresponding azimuth angle and zenith angle can be calculated by a slight adjustment of the proposed method ( Figure 4). The coordinates of the lunar surface are directly given in the PA. The following analyses use these criteria, unless otherwise specified.
Then, the coordinates of the intersection point and corresponding azimuth angle and zenith angle can be calculated by a slight adjustment of the proposed method ( Figure 4). The coordinates of the lunar surface are directly given in the PA. The following analyses use these criteria, unless otherwise specified.
Then, the coordinates of the intersection point and corresponding azimuth angle and zenith angle can be calculated by a slight adjustment of the proposed method ( Figure 4). The coordinates of the lunar surface are directly given in the PA. The following analyses use these criteria, unless otherwise specified. A qualitative comparison is shown in Table 1. The time is given in UTC +8. The STK results are the views from lunar center, while the images of Chang'E-3 and the simulation results of this paper are viewing from the lander. It can be inferred that the twilight arcs from three sources are largely the same (note that the attitude of the Earth and the rendering effect are slightly different). The disk views from the STK and this paper are similar, but the disk views from the Chang'E-3 lander are unidentifiable due to the low spatial resolution and high cloud coverage. Figure 5 shows infrared global geostationary composite v1 at UTC 2013/12/24 18:15. The cloud texture in the red crescent is very close to that in the Chang'E image acquired four minutes later.

13
MER.DE421 [38], i.e., (44.1412, −19.5382) in PA.DE430. The orientation angles differ slightly from the angles in Equation (7) [32]. In addition, the simulation results from the STK 11 are introduced as mutual authentication data. A more credible validation method will only rely on field measurements in future missions.
A qualitative comparison is shown in Table 1. The time is given in UTC +8. The STK results are the views from lunar center, while the images of Chang'E-3 and the simulation results of this paper are viewing from the lander. It can be inferred that the twilight arcs from three sources are largely the same (note that the attitude of the Earth and the rendering effect are slightly different). The disk views from the STK and this paper are similar, but the disk views from the Chang'E-3 lander are unidentifiable due to the low spatial resolution and high cloud coverage. Figure 5 shows infrared global geostationary composite v1 at UTC 2013/12/24 18:15.
The cloud texture in the red crescent is very close to that in the Chang'E image acquired four minutes later.

A quantitative comparison between the intermediate result and the STK outputs is also made.
Five groups of random numbers corresponding to the timestamp, lunar latitude, lunar longitude, azimuth bias, and zenith bias are combined to generate 100 check points. The aberration correction and light time compensation are not made in this comparison to keep pace with the software. Because the STK does not give built-in functions to report the results for a Moon-based sensor, we estimate the intersection position manually in the 3D interface point by point by using a conic beam as narrow as 2 × 10 −9 degrees. This method is very time-consuming. The comparison results are shown in Figure 6. All geo-distance deviations are smaller than 400 m and have no significant correlation to any parameters. The mean value and standard deviation are 108.6 m and 95.6 m, respectively. Interestingly, the data between 2009/01/01 and 2016/01/01 obviously fit much better than the data from other times (between two blue dashed lines). Even from 2010/01/01 to 2015/07/01, all the differences are smaller than 2 m. We think that this result may be related to the release date of the used software version and the operation of the ICRF since 2009/01/01, according to the STK help document. 14 smaller than 400 m and have no significant correlation to any parameters. The mean value and standard deviation are 108.6 m and 95.6 m, respectively. Interestingly, the data between 2009/01/01 and 2016/01/01 obviously fit much better than the data from other times (between two blue dashed lines). Even from 2010/01/01 to 2015/07/01, all the differences are smaller than 2 m. We think that this result may be related to the release date of the used software version and the operation of the ICRF since 2009/01/01, according to the STK help document.

The Scale of the Aberration Effect
t, and (4) with both treatments. The sensor is located at (−60,0). The differences between cases (1-3) and case (4) reveal the scale of these two effects. In Figure 7

The Scale of the Aberration Effect
According to Equation (18), the typical deviation of a moon-based sensor is about 1.3 km on the Earth's surface. The light time in Equation (23) is approximately 1.27 s and therefore leads to a 0.6 km deviation in geolocation. In most cases, as a result of the homodromy of the Earth's rotation and the Moon's revolution, the deviations cancel each other out, leaving a 0.7 km remainder. More detailed results are given in Figure 7. We first use Equation (27) to calculate the zenith angle and azimuth angle on the basis of the flow chart in Figure 4 and then reverse the procedure to recalculate the LOS and intersection points (1) when the aberration is not corrected, (2) when the light time is not compensated, (3) without either treatment, and (4) with both treatments. The sensor is located at (−60, 0). The differences between cases (1-3) and case (4) reveal the scale of these two effects. In Figure 7, the red dotted line represents the aberration effect, the blue dashed line represents the light time effect, and the black solid line represents the comprehensive effect. The geo-distance between corresponding points, used as an index, is in line with the analysis in the initial part of this paragraph. It should be noted that the comprehensive effect is weakly related to the seleno-graphic coordinates of the sensor but strongly connected with the position of the reflector. For example, if the reflector is near the edge of the Earth disk, the comprehensive effect can be larger than 10 km. One reason for this is that the same aberration angle will result in a larger projected geo-distance at the edge of Earth disk than at the disk center; in addition, the light time only causes small geo-distance deviations at high latitudes.

Pointing Deviation
There are some parameter uncertainties. We describe the general relationships between the errors of different parameters and their corresponding geo-distance deviations with examples from 2017 and the maximum and minimum deviations in the past 50 years. The LOS is pointed to the geocenter, and the relationship between the parameter errors and geo-distance deviations is approximately linear around the Earth's disk center. The aberration effect and light time effect of Moon-based Earth observation are two orders of magnitude larger than those of satellite imaging. Although the pointing accuracy requirements of passive optic sensors are not very high, they are considerably higher than the model accuracy analyzed in Section 5.3. Moreover, the compensation is essential if the model is used as a direct geo-location or correction model for array spectrometry image.

Pointing Deviation
There are some parameter uncertainties. We describe the general relationships between the errors of different parameters and their corresponding geo-distance deviations with examples from 2017 and the maximum and minimum deviations in the past 50 years. The LOS is pointed to the geocenter, and the relationship between the parameter errors and geo-distance deviations is approximately linear around the Earth's disk center.
Roughly, one second angular error will cause a geo-distance deviation of 1.84 km if the distance is 380,000 km. The spatial error will be directly projected on the Earth's surface with a scaling factor referring to the projection angle and Earth's curvature. The pointing errors attributed to some sources over time are shown in Figure 8. The figure only shows results from 2017, when the sensor was at (−60, 0); however, the geo-distance deviations of the other five sites over the past 50 years have also been investigated (Table 2). Generally, an angular error of one second is likely to generate a geo-distance deviation less than 1.94 km, while a 10 m positional error is likely to induce a geo-distance deviation less than 10 m. The ground deviation caused by the errors in zenith angle or lunar longitude are relatively not sensitive to time and seleno-graphic coordinate. They have more stable impact than the errors in other designable parameters. The positional errors can be converted to equivalent angular errors. The conversion refers to the nearest Earth-Moon distance to make sure that the maximum value is given. Therefore, the model accuracy can be quantitatively restricted by the sum of all kinds of angular errors.

Pointing Deviation
There are some parameter uncertainties. We describe the general relationships between the errors of different parameters and their corresponding geo-distance deviations with examples from 2017 and the maximum and minimum deviations in the past 50 years. The LOS is pointed to the geocenter, and the relationship between the parameter errors and geo-distance deviations is approximately linear around the Earth's disk center.  Table 2. The max and min geo-distance deviation (meters) derived from one second sensor attitude or seleno-graphic coordinate errors over the past 50 years when the line of sight (LOS) is pointed to the geocenter. The simulation step size is one day. The deviations caused by the librations and lunar center coordinates given by DE430 are consistent with Figure 8 and are therefore not shown in this .5-10 Roughly, one second angular error will cause a geo-distance deviation of 1.84 km if the distance is 380000 km. The spatial error will be directly projected on the Earth's surface with a scaling factor referring to the projection angle and Earth's curvature. The pointing errors attributed to some sources over time are shown in Figure 8. The figure only shows results from 2017, when the sensor was at (−60,0); however, the geo-distance deviations of the other five sites over the past 50 years have also been investigated (Table 2). Generally, an angular error of one second is likely to generate a geo-distance deviation less than 1.94 km, while a 10 m positional error is likely to induce a geo-distance deviation less than 10 m. The ground deviation caused by the errors in zenith angle or lunar longitude are relatively not sensitive to time and seleno-graphic coordinate. They have more stable impact than the errors in other  Table 2. The max and min geo-distance deviation (meters) derived from one second sensor attitude or seleno-graphic coordinate errors over the past 50 years when the line of sight (LOS) is pointed to the geocenter. The simulation step size is one day. The deviations caused by the librations and lunar center coordinates given by DE430 are consistent with Figure 8 and are therefore not shown in this

Model Accuracy
The error sources are divided into seven categories: Orbit, libration, seleno-graphic, pointing, propagation, EOP, and geographic ( Figure 9). The maximum error of each source will be evaluated.
In DE430, libration errors are coupled with orbit errors of approximately 0.1 . For the DE430 gravity field and librations, the uncertainty of the retroreflector PA coordinates is only 0.12-0.27 m [24]. If the error is 0.27 m, the deviation angle is approximately 0.032 . The coordinate accuracy of other sites is much lower. The coordinates of the Apollo lunar surface experiment packages may have errors as large as 30 m [39]. For cartography, the widely used Unified Lunar Control Network 2005 (ULCN 2005) is based on Clementine images acquired in 1994 and the previous ULCN. The solution root mean square is 0.9 pixels, corresponding to 100 m-300 m on the lunar surface [40]. The accuracy of next generation networks may reach 20 m at the resultant tie points [41]. For topography, the Lunar Reconnaissance Orbiter (LRO) Lunar Orbiter Laser Altimeter (LOLA) Digital Elevation Model (DEM) represents more than 6.5 billion measurements, and the average accuracy of each point after crossover correction is better than 20 m horizontally and approximately 1 m radially [42]. Therefore, the lunar map coordinate accuracy can be considered 20 m (2.4 ) for our simulation studies. After settlement on the Moon, the observatory position accuracy is expected to reach 0.2 m. 17 designable parameters. The positional errors can be converted to equivalent angular errors. The conversion refers to the nearest Earth-Moon distance to make sure that the maximum value is given. Therefore, the model accuracy can be quantitatively restricted by the sum of all kinds of angular errors.

Model Accuracy
The error sources are divided into seven categories: Orbit, libration, seleno-graphic, pointing, propagation, EOP, and geographic ( Figure 9). The maximum error of each source will be evaluated. m-300 m on the lunar surface [40]. The accuracy of next generation networks may reach 20 m at the resultant tie points [41]. For topography, the Lunar Reconnaissance Orbiter (LRO) Lunar Orbiter Laser Altimeter (LOLA) Digital Elevation Model (DEM) represents more than 6.5 billion measurements, and the average accuracy of each point after crossover correction is better than 20 m horizontally and approximately 1 m radially [42]. Therefore, the lunar map coordinate accuracy can be considered 20 m (2.4″) for our simulation studies. After settlement on the Moon, the observatory position accuracy is expected to reach 0.2 m.
In astronomy, a telescope or radio telescope is required to point to a star or radio source whose right ascension and declination position is already known. The position of a star or radio source seldom changes, so errors are mainly derived from the pointing accuracy of the telescope. At present, the pointing accuracy of a telescope after correction is typically 1″on the Earth's surface [31,43]. Without atmospheric distortion and with slight crustal deformation variation, the pointing accuracy of a Moon-based sensor may be higher.
The atmospheric refraction is related to wavelength. For visible light, the ground deviation In astronomy, a telescope or radio telescope is required to point to a star or radio source whose right ascension and declination position is already known. The position of a star or radio source seldom changes, so errors are mainly derived from the pointing accuracy of the telescope. At present, the pointing accuracy of a telescope after correction is typically 1 on the Earth's surface [31,43]. Without atmospheric distortion and with slight crustal deformation variation, the pointing accuracy of a Moon-based sensor may be higher.
The atmospheric refraction is related to wavelength. For visible light, the ground deviation is usually 5 m (0.003 ) when the local incident angle is 45 • [34]. This deviation will be larger at a large incident angle. For microwaves, the ionospheric effect is more significant. The ground deviation can reach 100 m in the ultrahigh frequency band [44]. The stratospheric effects are much smaller. Thus, we assume the total refractive deviation is 150 m, corresponding to 0.086 .
The EOP and geographic errors usually refer to the geo-center, and should be converted to refer to the moon-based sensors in our analysis. In Moon-based Earth observation, the target position in the right ascension and declination relative to the sensor is always changing, so the transformations introduce extra errors. The model should be reversed to calculate the right ascension and declination target position. Rewrite Equation (23) as: where t 1 and is already given. An iteration is used to solve the equation. The initial error is 10 m or 0.006 , and almost no estimated error remains after several iterations. However, the SOFA algorithm to calculate the orientation matrix has a max error of approximately 0.9 in 21st century [45]. Relative to a Moon-based sensor, this error is less than 0.016 .
The impact of the Earth's DEM can be easily corrected. In the forward model, the intersection with the DEM should be numerically coupled with Equation (24). In the reversed model, the elevation is already contained in −−−−−−−→ A ITRS (t 1 ). However, each DEM contains intrinsic errors in relation with a particular terrain and land cover type [46]. The global assessment result shows that for Shuttle Radar Topography Mission (SRTM) DEM the absolute geolocation error is less than 15 m and the absolute height error is less than 10 m [47]. Thus, the geographic error is expected less than 25 m, or 0.015 referring to the sensors.
The estimated accuracies are summarized in Table 3. All other errors are assumed to be within 0.1 . The main errors are the pointing control error and the position error on lunar surface. A great boost on seleno-graphic accuracy is expectable after the settlement on the Moon. At that time, the lunar surface positioning error will be as small as the position error of the LLR reflectors, and the overall model accuracy will be better than 1.5 . It is high enough for mission operation and image geolocation. At this current stage, if we use the LRO LOLA DEM to simulate data, the overall model accuracy is about 4 . Even in this case, the accuracy is acceptable compared with the limitations presented at the end of Section 2.

Predicting the Pointing Direction and the Observational Constraints of Site Selection
In passive optic sensor cases, to point the LOS to the geocenter, the sensor attitude is determined by only the location of the observatory and the height of the sensor. Figure 10 shows the change rule of the azimuth angle and zenith angle when the height is set to 0. The seleno-graphic coordinates of the sensors in subfigures (a), (b), and (c) are (0,0), (0,80), and (80,0), respectively. The black solid line shows the zenith angle, the red dotted line shows the azimuth angle and the blue dashed line shows the angle between the LOS and the negative light direction toward the sensor (LOS-SUN angle). The negative zenith angle results from the definition of the rotation direction. As stated in Section 5.3, the overall prediction accuracy is expected to be less than 1.5 for mission operation.
The LOS-SUN curve in Figure 10 indicates the risk of the sunlight going directly into the sensor. The closer the curve is to 0, the higher the risk is. The similarity of the curves in different subfigures suggests that the LOS-SUN angle is only slightly affected by the sensor's position. Thus, its monthly approach to 0 is inevitable. To reduce the probability of the Sun entering the FOV, we suggest the use of a small FOV that is slightly larger than the Earth disk.
The irregular relative motion causes the projection of the Earth on the lunar surface to move around the seleno-graphic origin, and the selenocentric angle is typically less than 10 • . If the sensor is put in this area (e.g., a circle around the seleno-graphic origin with a 10 • radius), the azimuth angle must vary from 0 • to 360 • in each month to track the Earth (Figure 10a). Otherwise, the azimuth angle will be biased to a smaller range, and the farther the sensor is from the origin, the slower the angle will vary ( Figure 10). This necessitates a low-frequency operation of the turntable. In passive optic sensor cases, to point the LOS to the geocenter, the sensor attitude is determined by only the location of the observatory and the height of the sensor. Figure 10 shows the change rule of the azimuth angle and zenith angle when the height is set to 0.  Unlike the azimuth angle, no matter where the sensor is put, the zenith angle always varies within a range of a dozen degrees. However, if the absolute value of the zenith angle is too large, the signal from the Earth will be contaminated by the reflected or emitted energy from the lunar surface. To avoid this contamination, the absolute value of the zenith angle should be smaller than 90 • − FOV/2. In Figure 10b,c, the absolute value of the zenith angle is larger than 85 • in some time periods. Therefore, places outside the 80 • circle should be carefully examined. Figure 11 shows a 19-year simulation on how the observatory position impacts the observational parameters, including the percentage of full coverage, the mean derivative of elevation angle and azimuth angle, and the number of times when LOS-SUN angle is smaller than FOV/2. As many former studies shows, the full coverage area, where the FOV are not sheltered all the time, is almost in an 81-degree circle around the seleno-graphic origin on lunar surface. At the edge of the Moon's nearside, the full coverage rate decreases to 0.4 (Figure 11a). That means the sensors should be put inside the 81-degree circle to access to full coverage data at any time. Figure 11b,c shows the derivative of the two pointing angles. The change of elevation angle is slow but the azimuth changes intensively near the seleno-graphic center. The azimuth adjustment rate will increase 10 times if the observatory is moved from the edge of the Moon's nearside to the center. That will result an extra burden for SAR pointing system. To minimize the pointing adjustment, the observatory should be put far from the center, e.g., areas close to the 81 • circle on the lunar near side hemisphere. Compared with the south pole regions, the north polar regions are more likely to be affected by the direct sunlight. The best places to avoid sunlight (the red circle in Figure 11d) is the north tropical zone on the west part of the lunar nearside, about 5% better than the north polar regions. Therefore, after the comprehensive analysis of four subfigures, (15, −75), near the Crater Krafft, may be one of the best sites for our moon-based earth observatory, due to a relatively low possibility of direct sunlit, slow adjustment of pointing direction, and full coverage of the earth disk all through the time. The south pole area is considered a hot candidate for the permanent Moon base [48,49].
Because of their high altitude, some regions near the south pole have good Earth visibility [50].
On the other hand, since the outer rim of the South Pole-Aitken basin at the near side can be seen from the Earth, this region may provide a suitable place for the attached observatory [51].
According to Figure 11, the area near (-81,0) are not the best but also a quite good place. Another The south pole area is considered a hot candidate for the permanent Moon base [48,49]. Because of their high altitude, some regions near the south pole have good Earth visibility [50]. On the other hand, since the outer rim of the South Pole-Aitken basin at the near side can be seen from the Earth, this region may provide a suitable place for the attached observatory [51]. According to Figure 11, the area near (−81,0) are not the best but also a quite good place. Another crucial factor in this area for the visibility analysis is the complex terrain. More detailed analyses require high-precision local DEM, and thus will be reported in following papers.

Simulating the Range History of Moon-Base SAR
The SAR imaging mechanism is quite different from passive optic sensors, but the relative motion patterns are the same. Therefore, we can build the moon-based SAR geometry from modifying the proposed geometric model.
For the sake of light time effect, the 'stop-go' hypothesis is no longer valid [18,52]. The SAR baseband signal of a point target echo after demodulation can be written as [39,53]: In this equation, τ denotes the fast time taking the center of each pulse as the origin. t denotes the slow time referring to UTC. A 0 is a complex constant related to the properties of the point target. The range window ω r refers to the pulse envelop. The one-way azimuth window p a (θ) ≈ sin c (θL a f 0 /c), where L a is the antenna length. ϑ(t) is the angle between the wave and the normal plain of the antenna in range direction. f 0 is the carrier frequency and K r is the chirp rate of the transmitting signal. To simulate the signal, the propagation range history R p (t) must be accurately determined. Figure 12 shows the round-trip propagation geometry of an orbit SAR system. In a SAR system, the transmitting time t can be recorded, while the 'go and back' propagation time ∆t 1 and ∆t 2 is unknow. Considering the light time effects, we assume that the pulse center leaves the antenna S at t, arrives at the target A at t + ∆t 1 , and reaches the antenna S again at t + ∆t 1 + ∆t 2 . Then the propagation range is: where An iteration method or an approximation similar to Section 3.7 can be used to numerically solve Equation (31). The numerical result from Equations (29)- (31) gives an accurate signal simulation which can be used to develop fast focusing methods.

S(t+∆t 1 +∆t 2 ) A(t)
A(t+∆t 1  In SAR footprint geolocation, we need two targets to determine the three pointing angles. The first constrains the LOS, and the second describes the third pointing angle in Equation (5). For example, if the long side of the antenna is required to approximately parallel to the Earth's longitude, the constraint condition can be written as: where A and B indicate two targets with same longitude; x, y, and α R in M SCS_OLCS are the three unknown numbers.
The aberration effect will result in a small offset angle in both range and azimuth. Theoretically, the transmitting wave and receiving wave should be translated to SCS to calculate the transmitting and receiving pattern weights. However, concerning the typical aberration angle is only 3.4 µrad, the difference is relatively small. For example, assuming the wavelength is 0.03 m (X-band) and the antenna length is 60 m, the round-trip weight difference is less than 3% in the half power beam width.
The angular aberration effect almost gives no impact on range history, but the light time will cause slight changes in both transmitting and receiving propagation range. The real range history should be the sum of (I) transmitting  [18]. In airborne SAR and low earth orbit SAR imaging, the range history is approximated to (III) + (III). In reference [18], it is approximated to (III) + (IV). The differences between (III) + (III) and (I) + (II) reach 44.41 m, 34.12 m, and 7.90 m within a 10-min synthetic aperture when the target is put from low latitude to high latitude places ( Figure 13). The variation of difference between (III) + (IV) and (I) + (II) are 22.21 m, 17.06 m, and 3.95 m. Though the residual error reduced by half, it is still greater than the signal wavelength during the long synthetic aperture time (SAT). The detrended results suggest that the variation of residual errors can be described as a linear term (Figure 13d). The high order error is within 1 mm in 200 s SAT and 3 mm in 600s SAT around the time of closest range, at least one order smaller than the wavelength of L, C, and X band. That means though the approximation (III) + (III) or (III) + (IV) is not accurate enough, obvious defocusing will not arise. The linear errors can only cause a mismatch on focusing position.

23
The aberration effect will result in a small offset angle in both range and azimuth. Theoretically, the transmitting wave and receiving wave should be translated to SCS to calculate the transmitting and receiving pattern weights. However, concerning the typical aberration angle is only 3.4 μrad, the difference is relatively small. For example, assuming the wavelength is 0.03 m (X-band) and the antenna length is 60 m, the round-trip weight difference is less than 3% in the half power beam width. The angular aberration effect almost gives no impact on range history, but the light time will cause slight changes in both transmitting and receiving propagation range. The real range

Conclusions
The correspondence between the intersection point and sensor attitude is given by the proposed geometric model. The simulated disk view is in good accordance with the image from the Chang'E lander and STK output. The analysis of the forecasted sensor's attitude suggests that most places at the near side of the Moon are geometrically suitable for Earth observation, except for the marginal areas of the lunar disk and the neighborhood around the seleno-graphic origin. An observatory at the edge of the disk will likely experience view blocking by the topography or even by the lunar sphere itself, while an observatory near the lunar disk center will have a relatively high frequency of operation. The periodical approach between the Earth and the Sun on the celestial sphere requires a small FOV to prevent direct radiation. The analyses show that the areas near Cater Krafft have relatively low possibility of direct sunlit.
The light aberration scale is approximately 3.4 × 10 −6 rad, corresponding to a deviation of 1.3 km on the Earth's surface, while the light time scale is approximately 0.6 km in the opposite direction. As a result of neutralization, the combined effect is approximately 0.7 km, which is in the same order of magnitude as the pointing error. The light time also causes an error greater that the wavelength in moon-based SAR range history. In a 600s SAT around the time of closest range, the linear error can be as large as 44 m, but the high order error is at least one order smaller than the wavelength. The overall model accuracy is evaluated to be better than 4 at current stage and will increase to be better than 1.5 when the lunar positioning error is restricted to the level of the LLR reflector seleno-location accuracy.