CLUSIM: A Synthetic Aperture Radar Clutter Simulator for Planetary Exploration

In this paper we present the CLUtter SIMulator (CLUSIM), a special program simulating radar side echoes from rough planetary surfaces using realistic topography data sets. A numerical model of realistic topography of the Martian surface, based on Mars Orbiter Laser Altimeter data, is developed. A specially developed computer routine for evaluation of wideband radar echoes reflected from rough surfaces, capable of aperture synthesis simulation, is described. A synthetic radargram for a portion of Mars Express (MEX) orbit 9466 is computed and validated against experimental data obtained by the MARSIS radar instrument. Finally, a previously developed ionospheric phase correction procedure is numerically tested with new simulated echo signals. Impact of the surface clutter on the ionospheric correction procedure is investigated with a direct numerical comparison to a known benchmark result, which shows robustness of the correction algorithm with respect to the surface clutter.

A monopole antenna, with zero gain in the nadir direction, receives side clutter only. This was then subtracted from the main dipole antenna signal. However, complete elimination of the clutter echoes is not yet possible. More sophisticated radar techniques, such as focused aperture synthesis and synthetic aperture radar (SAR) interferometry, which require much higher data transfer and processing rates, might be a solution of this problem in the future. Now surface clutter can be simulated, compared to real experiment data, and subtracted from it more or less accurately. Thus, numerical simulation of the clutter is a necessary stage of the radar data analysis and interpretation procedure (Nunes et al., 2010(Nunes et al., , 2011. Theoretical clutter estimates are also required at all the stages of the radar experiment, including instrument design and observational strategy planning. For these purposes, one needs a radar experiment model, allowing for these estimates. This model in turn should consist of the physical model of the object structure (surface) and algorithm for the numerical solution of the electromagnetic equations, ruling an interaction of the electromagnetic wave with the object.
Development of SAR clutter simulation tool, capable of high-resolution topographic data assimilation and efficiently exploiting powerful resources of modern computing machinery available now for the Martian applications, is the motivation of this investigation.
The paper is organized as follows. The new clutter simulator CLUSIM, which is the main goal of this study, is described in section 2 (the Martian surface topography model and the algorithm for the electromagnetic scattering simulation in sections 2.1 and 2.2, respectively). Section 3 is entirely devoted to the ionospheric correction algorithm test, including the Martian ionospheric model and phase correction approach. Simulation results are validated against the experimentally measured data from the MARSIS instrument and briefly discussed in section 5, and in section 6 conlusions are formulated and some final remarks are given.

Synthetic Aperture Radar Clutter Simulator
Clutter simulators (Nouvel et al., 2004;Spagnuolo et al., 2011;Russo et al., 2008) consist of two basic structure units: a topography interpolator and a radar echo simulator itself.
Both the rough surface model and the approach to the electromagnetic equations itself can be deterministic or statistical. To validate the detection of subsurface interfaces, numerical electromagnetic models of surface scattering should be used to produce simulations of surface echoes, which are then compared to real echoes detected by the radar. A realistic model of the surface should be based on real observations. Accounting for these considerations, a special package of computer codes named CLUSIM (CLUtter SIMulator) has been developed and tested by us for the Mars case. The package consists of the topography analyzer and numerical integration routine. Topography analyzer block in turn includes routines for Mars Orbiter Laser Altimeter (MOLA) data assimilation, their preprocessing and preselection, removal of suspicious (nonreliable) data, topography interpolation, and visual control.
The main purpose of the package is the simulation of radargrams at any given part of the spacecraft orbit.
For Mars, the MOLA database is probably the most used and best developed source of surface topography data. These data are available in the form of structured list of latitudes, longitudes, and heights of the MOLA spots grouped in orbital tracks. The spatial resolution of these data is about 300 m along a track and up to several kilometers across it from track to track (Smith et al., 2001). To build a digital elevation model (DEM), which is an important part of the clutter simulating algorithm, these data must be interpolated on the planetary surface, giving the surface height as a function of latitude and longitude.
To do this, the following approach is applied. All the available MOLA spots are projected on the tangent plane to the subsatellite point of the planetary surface, keeping their ordered track structure. After that, the height profile along each track is approximated by a cubic spline. Then, a regular grid of parallel lines is built on the tangent surface plane. This grid also considers the grid lines of the height profile, which are defined by the MOLA tracks in the crossing points. Again, this is approximated by a cubic spline. Finally, a uniformly spaced rectangular grid is built on these grid lines, which gives an approximation of the rough surface terrain. Accuracy of this approximation in a given point of the surface should be on the order of height difference between adjacent MOLA spots nearest to this given point. On the MOLA spots, this approximation exactly corresponds to the MOLA surface elevation data. Thus, the surface topography is defined as a function of Cartesian coordinates on the tangent plane. The curvature of the spherical surface (deviations of the spherical surface from the tangent plane) can be now just added to the topography approximation, if necessary. In practice, the impact of the planetary surface curvature is small, so the tangent plane approximations yield the satisfactory results.
With that approximate topography defined on the tangent plane with Cartesian coordinates, one can evaluate the numerical solution of the electromagnetic problem of interaction of the sounding wave with the object of known shape.

Algorithms for Solving Electromagnetic Field Equations
An approach to the numerical solving of the electromagnetic waves equations can also be deterministic (Fa & Jin, 2010;Liu et al., 2014) or statistical (Ilyushin, 2014). In the deterministic case of an object of a given structure, there are no principal difficulties to numerically solve electrodynamical equations. A number of algorithms, such as finite difference time domain (FDTD) (Yee, 1966), finite element method (FEM) (Zienkiewicz et al., 2013), discrete dipole approximation (DDA) (Purcell & Pennypacker, 1973), and T-matrix (Waterman, 1965) are well developed now.
However, planetary radar echo simulations are characterized by a large object size (compared to the wavelength) and by a wide frequency spectrum of the signal. Strict approaches of the computational electrodynamics mentioned above (FEM, FDTD, DDA, etc.), are in fact inapplicable to this particular case of the SAR sounding of the planetary surface because of the large size of the object. The known reported FDTD simulations are one dimensional or two dimensional (Xu et al., 2006) or restricted to a little domain of the medium (Boisson et al., 2009). Application of the aperture synthesis and other averaging signal processing techniques, such as migration (Herique & Kofman, 1997), requires their adequate representation in the numerical model, which in turn results in many repeated calculations of the same type with different input parameters, and the computation load is vastly increased. Estimates of mean signal parameters, averaged over many realizations, result in additional increase of computational complexity. Probably, for this reason, most known clutter simulations are more or less simplified, for example, they lack the aperture synthesis (Fa & Jin, 2010;Smirnov et al., 2014). In other simulation algorithms, the aperture synthesis is performed at the postprocessing stages (Berquin et al., 2015).
Thus, for such problems of scattering from the large rough surfaces, approximate techniques are more appropriate. Since the height profile is not small compared to the wavelength, the small perturbation technique is also not appropriate. The only feasible approach is Kirchhoff approximation (Beckmann & Spizzichino, 1963;Ogilvy, 1991), which is typically used to solve these kinds of problem. Making use of the vectorial Huygens-Fresnel principle, which is a counterpart for the conventional scalar Kirchoff approximation, allows one to consider polarization of the radar echoes and depolarizing effects. However, depolarization is known to be absent within the gently undulating surface approximation (Beckmann & Spizzichino, 1963;Ogilvy, 1991), which is a reasonable assumption for many planetary terrains. On the other hand, much depolarization of the reflected waves comes from the discrete objects and small-scale structures, which available DEM of planetary surfaces do not account for.
It is worth mentioning here that given the wavelength of the MARSIS carrier signals, the presence of rocks, even if they were several decimeters in length, is unlikely to significantly affect the echo signal. Neumann et al. (2003) estimated the roughness of the Martian surface within the MOLA footprint (75 m across) from the widening of the pulse echo. They found that over the vast majority of the Martian surface, the root-mean-square (RMS) height of the surface is at most a few meters. Notable exceptions are the Olympus Mons areole and the Olympia Undae dune field. However, for much of the Martian surface, the small-scale roughness height is a very small fraction of the MARSIS wavelength, which ranges between 60 and 170 m.

10.1002/2017RS006265
In addition, a single reflection from a rough surface-to-vacuum interface, typical for the terrestrial planets, does not produce high reflection amplitude nor significant change of the wave polarization state in the experimental planetary radar studies (Black, Campbell, & Ostro, 2001;Hagfors, 1967). Both anomalously high radar albedo and strong depolarization of radar pulses are typical for Galilean satellites of Jupiter (Ostro & Shoemaker, 1990). These remarkable properties of the Jovian icy moons have been attributed to the volume scattering (Hagfors et al., 1997), multiple subsurface reflections (Eshleman, 1986;Ostro & Pettengill, 1978), and coherent backscattering enhancement effect (Black, Campbell, & Nicholson, 2001;Hapke, 1990;Hapke & Blewett, 1991). Since both volume scattering and coherent backscattering enhancement are not expected in the Martian crust at MARSIS wavelengths, we do not include these effects in our model. These effects are discussed in Ilyushin (2012). In our present study, we neglect the polarization of the waves and restrict our consideration to the conventional scalar Kirchhoff approximation. As it has been shown by Ilyushin (2014), the calculation of the reflected radar signal spectrum (complex surface reflection coefficient R( )) within this approximation reduces to the numerical evaluation of the Kirchhoff integral over the whole reflecting surface: where r 0(t) and r are the source (spacecraft) location point and surface point of integration, respectively, and h(x, y) is the local surface height.
The Green function of the Helmholtz equation in the paraxial approximation is equal to where k = ∕c (Goodman, 1996).
For a dielectric surface with finite permittivity, this expression should be multiplied by the Fresnel reflection coefficient R F . Due to the near-nadir sounding geometry of the orbital GPR instruments and moderate slopes typical for much of the Martian surface, the Fresnel reflection for normal incidence is a reasonable approximation where is a permittivity of the surface. Here we assume that the surface is homogeneous and its permittivity is uniform everywhere.
In the Kirchhoff approximation, the complex amplitude of the radar echo depends on the dielectric permittivity only through this coefficient (3). Since the primary objective of our simulations is the ratio of the main nadir echo to the side clutter echoes rather than signal-to-noise ratio (SNR), we do not pursue the absolute signal level calibration here and simply assume R F = 1 everywhere, which corresponds to an ideal reflective surface.
If the aperture synthesis is applied, an extra integration over the spacecraft coordinate is needed (see Figure 1). Aperture synthesis is essentially averaging over a portion of a satellite trajectory of length 2L (synthetic aperture) with a weighting function exp , where s is the coordinate along the trajectory, is the Doppler filter number (positive or negative integer). In this paper we consider only the most important case = 0.
Averaging the equation (1) over the synthetic aperture, we obtain assuming that satellite trajectory r 0(t) = x 0 (t), 0, z is parallel to the x axis and the horizontal. Small deviations from a horizontal trajectory can be easily accounted for in the phase exponent in (4). Integrals of the quadratic exponent involved here, resulting in the error function of complex argument in equation (4), can be effectively expressed through the so-called Faddeeva function (Abramowitz & Stegun, 1972). Recently, a computer code package in C++ and Fortran languages has been released by Johnson (2012). With it, all the computations can be implemented on modern high-performance computing machines. For the piecewise-planar approximation of the surface elevation with triangular of rectangular facets, the complex reflectance of an individual facet R f ( ) can be derived from (4) as follows (see Figure 2): where x 1 (y) and x 2 (y) are the left and right boundaries of the given facet. The inner integration (over the x variable) can be performed analytically over the triangular of rectangular facet (see Ilyushin, 2004). Total echo signal from the whole rough surface is therefore a sum of contributions of all the facets R f ( ), all of which is an one-fold integral in spatial domain. Berquin et al. (2015) make use of the vectorial Huygens-Fresnel principle and derive an expression for polarized radar echo from the planetary surface: with the dyadic Green functionḠ where E and H denote electric and magnetic fields, respectively,n is a unit normal to the surface, x is a point on the surface, = √ ∕ is the characteristic impedance of the surface, and are the permittivity and permeability, k is the wave number andk s is the unit vector pointing in the scattering direction. Assuming the length of the synthetic aperture to be small compared to the spacecraft orbit height, we can regardk s as a constant and integrate over the synthetic aperture analytically, as is explained above for the scalar case. However, if the gently undulating surface approximation is adopted, there are no depolarization of the radar echoes (Beckmann & Spizzichino, 1963;Ogilvy, 1991). On scales well exceeding the sounding wavelengths, for which the Kirchhoff approximation and Huygens-Fresnel principle are valid, most Martian terrain can be regarded as smooth, gently undulating surfaces (Kreslavsky & Head III, 2000). The radar echo can be more or less depolarized by the surface features (rocks, boulders, etc.) comparable with the wavelength. However, because planetary GPR instruments have not yet had polarimetric capabilities, for purely energetic assessments we perform our numerical simulations within the conventional scalar Kirchhoff approximation. Berquin et al. (2015) also derive an expression for the diffraction integral over the whole surface of the triangular facet as an approximate sum of an infinite series. These approximate expressions cannot be further integrated analytically, so the aperture synthesis can only be performed numerically on the postprocessing stage. In addition to complicated procedure of adding up many terms of the sum, the approximation error results in spurious reflections coming from the edges of the facets, in turn causing the Bragg diffraction maxima produced even by the perfectly flat surface approximated by such facets. To remove this, Berquin et al. (2015) proposed to use irregular grids for the facet approximation. Our approach with analytical aperture synthesis is free of these Bragg diffractional artifacts and provides more flexibility in the choice of DEM, since it is not confined to the constant phase and linear phase approximations of the surface. The aperture synthesis approach, which is proposed here, can also be applied to the case of vector field simulation. The exponent in (7) is expanded in series in small variations of the spacecraft coordinate x 0 and integrated over, yielding the expression containing the error function like (4) as well as in the scalar case. After that, further integration over the rough surface can be performed using different approaches, thus avoiding emulation of the aperture synthesis at the postprocessing stage.
When the radar echo spectrum is evaluated, the signal can be obtained with standard signal processing technique (matched filtration). An expression for the compressed radar signal s(t) with partial correction of the ionospheric phase distortions (Ilyushin & Kunitsyn, 2004) writes Input data (Martian topography from MOLA server) where = 2 f is the angular frequency, k = c is the wave number, H( ) is the spectral window function, F( ) is the initial radar pulse spectrum radiated by an instrument, R( ) is a complex reflection coefficient of rough surface, n(z) is the height profile of the ionospheric refractive index, p (z) is the angular plasma frequency, N(z) is the ionospheric plasma electron concentration, r e is the classical electron radius, and c is the light speed in vacuum. Natural ionospheric phase shift ( ) is adaptively compensated by the phase curvẽ( ) (Ilyushin & Kunitsyn, 2004). As a compensating phase term, a phase curve ( ) of artificially constructed model plasma layer with negative sign is often used. Ilyushin and Kunitsyn (2004) analyzed a wide set of radio occultation profiles of Martian ionosphere available at that time and found that the triangular vertical profile of plasma density distribution is optimal for compensating phase distortions occurring in real Martian ionosphere.

Along-track interpolation
Further radar signal processing methods (migration, etc.) can be then applied to the simulated signals if necessary. As it has been shown by Ilyushin, (2010Ilyushin, ( , 2009a, radar wave scattering due to the plasma frequency fluctuations is relatively weak, except for specific ionospheric situations (high plasma frequency and intense fluctuations during strong solar events). Safaeinili et al. (2003) studied not only ionospheric phase distortion but also impact of the effect of Faraday rotation of wave polarization plane. They showed that in weak local magnetic fields, typical for Mars, this effect is not pronounced and may be neglected, even in measurements with polarization-sensitive antennas. Thus, in most situations, the ionosphere can be regarded as a flat layered slab of an isotropic dispersive medium. Under assumed approximations the ionospheric phase delay ( ) and the phase correction term̃( ) are represented by separate factors and can thus be considered independently in the simulations. We omit it everywhere, except in the section describing a test of ionospheric correction algorithm with newly simulated clutter data. (parallel and perpendicular to the spacecraft trajectory), integration along this trajectory (aperture synthesis) can also be performed analytically, yielding a result expressed through the Fresnel integrals in closed form.
In some specific cases, practical realization of facet approach is rather effective and physically adequate. Ilyushin (2004) adopted a piecewise-planar approximation of Martian polar trough profile and simulated MARSIS radar echoes coming from these troughs for different mutual orientation of the polar troughs and spacecraft trajectory. Essentially, the one-dimensional character of the valley allowed for continuous piecewise linear profile, ensuring high simulation accuracy at reasonable computational expenses. It has been shown by Ilyushin (2004) that for proper (normal) mutual orientation of the orbit and the trough, surface clutter generated by the trough is well suppressed by aperture synthesis (practically, below the dynamic range of the radar sounder).
However, the rectangular plane tile (facet) approximation of arbitrary 3-D rough surface is, generally speaking, discontinuous ( Figure 3). Such discontinuities of the model surface should generate artifacts, contaminating the result with unphysical reflections.
Continuous approximations are more complicated. The simplest one is the triangulation (planar approximation on the triangular grid, Figure 4). Surface triangulation has been used for clutter simulation in some radar studies (Fa & Jin, 2010;Liu et al., 2014); however, aperture synthesis technique was not applied there. To do this, one has to average the complex reflection amplitude (1) over some portion of the satellite trajectory. (1) over the triangular domain cannot be expressed in closed analytical form. However, integration over the spacecraft trajectory can be analytically performed. Thus, evaluation of the signal spectrum reduces to the numerical estimate of the Kirchhoff diffraction integral (1) over the whole surface with the integration kernel (4). The compressed radar pulse in time domain can then be evaluated by common matched filtration technique. We used the same pulse compression procedure as in our previous studies, where this procedure is described in detail by Ilyushin, (2004Ilyushin, ( , 2008Ilyushin, ( , 2009b. Systematic phase shift, introduced by the ionosphere, can be introduced separately as an independent exponential factor at the signal compression stage.

Diffraction integral in equation
The structure of the clutter simulation algorithm is shown in the flowchart ( Figure 5). Raw topographic data, obtained from the MOLA data archive (stage I), are interpolated along each track (stage II) and then across the tracks (stage III). For each point of the orbit, where the sounding with a sequence of chirp pulses (synthetic aperture) is performed, the interpolated topographic data are triangulated on the regular grid centered around the subsatellite point (stage IV). After that, the complex reflection coefficient of the rough surface (4) is evaluated (stage V) by the integration over the surface for each of the discrete frequencies, constituting the digital representation of the radar signal spectrum, in the cycle (stage VI). In the simulation 512 frequencies are used, as well as in the current MARSIS experiment.  After that, the signal is evaluated from its spectrum with the fast Fourier transform (FFT) procedure (stage VII), using common matched filtration procedure with Hann spectral window (Harris, 1978), according to the MARSIS technical specifications (Picardi et al., 1999;Calabrese et al., 2001). The signals are routinely calculated for a number of points along the MEX spacecraft orbit (stage VIII, cycle) for which the measured signal records are available. All the signals are then assembled into the radargram (stage IX, final), characterizing radar clutter echoes coming from the surface terrain on the given portion of the planetary surface landscape. The radargram is then represented graphically in the standard way.

Ionospheric Model and Correction Algorithm Tests
As planetary atmospheres are ionized by solar radiation, they become a dispersive medium for relatively long radar waves, used in subsurface sounding. The radar signals become distorted by the ionospheric phase dispersion (Cartacci et al., 2013;Safaeinili et al., 2003;Sanchez-Cano et al., 2015). The radar data processing routine must therefore include some special procedure for extra phase shift compensation (Armand & Smirnov, 2003;Ilyushin & Kunitsyn, 2004;Mouginot et al., 2008), adaptive or not. Impact of random wave scattering on the function of these algorithms have been previously studied (Ilyushin, 2008;Ilyushin et al., 2012). In this study, we perform a numerical test of the previously developed phase compensation algorithm (Ilyushin & Kunitsyn, 2004) together with a realistic surface clutter model. The reader is referred to Ilyushin and Kunitsyn (2004) for the detailed description of the ionospheric correction approach tested here.
For these simulations, we use a new developed empirical model for the dayside electron density of the Martian ionosphere (primary and secondary layers), called NeMars (Sanchez-Cano et al., 2013). The model is largely based on MARSIS Active Ionospheric Sounding (AIS) mode data (Gurnett et al., 2005;Morgan et al., 2013;Sánchez-Cano et al., 2012) and, to a lesser extent, on radio occultation data from the Mars Global Surveyor mission. In addition, this model partially assumes the Chapman theory. This model reproduces to a reasonable degree the main characteristics of the vertical electron density profiles obtained with the two techniques by considering solar zenith angle, solar flux F 10.7 as a proxy of the solar activity, and heliocentric distance. Typical model profiles used for the present simulation are shown in Figure 6.
For testing purposes, a special MARSIS radargram with adaptive compensation of model ionospheric distortions was simulated. As a source of initial ionospheric distortions, the NeMars model has been used. For the adaptive compensation of the distortion, a phase curve of a plasma slab layer with triangular height profile of the electron density (Ilyushin & Kunitsyn, 2004), has been used. Partial phase distortion is mitigated by adjusting the phase curve of the correcting plasma layer close to real ionospheric phase curve, which is unknown, and subtracting it from the signal phase. In practice, the parameters of correcting plasma layer (plasma critical frequency 2 c and ionospheric layer thickness H) have been adjusted for optimal (maximal) intensity contrast (Ilyushin & Kunitsyn, 2004) In this work, the adaptive ionospheric phase correction is performed on the post processing stage, after the matched filtration. However, within the paraxial approximation adopted by us, phase and clutter responses are multiplicative and can therefore be evaluated independently.

Results of Numerical Simulations
We present here sample results obtained for the southmost portion of MEX 9466 orbit. MOLA topography of the surface landscape surrounding the subsatellite track and its visual image provided by the High-Resolution Stereo Camera (HRSC) (Neukum & Jaumann, 2004) are shown in Figures 7  and 8, respectively. Figure 9 consists of five panels, each of them represents a measure or simulated radargram. Each radargram shows the intensity of radar echo in gray shades against echo signal delay (vertical axis) and position of the spacecraft in the orbit (horizontal axis) in the standard way typically used for radar data representation. Simulated radargrams are shown in Figures 9b-9e, together with the experimental MARSIS radargram (Figure 9a).
The radargrams presented in this paper represent about 118 s of observations at an altitude of 390 km and a speed of 4.2 km/s. MARSIS transmitted a linearly modulated, 350 s long pulse (a "chirp") that can be centered at 1.8, 3, 4, or 5 MHz, with a 1 MHz bandwidth. In the observations considered in this paper, the central frequency of the chirp was 5 MHz. The pulse echo is downconverted to a central frequency of 0.7 MHz, and sampled by the analog-to-digital converter at 2.8 MHz sampling rate. In-phase and quadrature components (I/Q) synthesis then reduces the sampling rate at 1.4 MHz. The synthetic aperture length 2L during observations ranged from 2,016 to 2,226 m. In the simulation, the value 2L = 2, 500 m was used, which is close to the optimal value for the unfocused aperture synthesis (half Fresnel zone size √ z∕2 = 2, 100 m).
For the assessment of the aperture synthesis efficiency, a single-pulse radar sounder echo has been simulated (shown in Figure 10 together with the synthetic aperture echo, for comparison). Namely, solid and dotted curves in Figure 10 represent radar signals with the application of the aperture synthesis and without it.
Finally, to comparatively investigate complexity of ionospheric compensation with ideally flat and realistic rough surface, we show the plots of intensity contrast C 2 I (13) of simulated radar signals with partial compensation of the ionospheric phase distortion. The graphs of the contrast calculated for the reflection from ideally flat surface and realistic Martian surface terrain are shown in Figures 11 and 12, respectively. There is shown dependence of the contrast function C 2 I on the two leading terms of the Taylor series of phase mismatch ( ) −̃( ) in (8), that is, a 2 ( − 0 ) 2 and a 3 ( − 0 ) 3 , where 0 is the central angular frequency of the radar chirp band.

Discussion of the Numerical Results
Numerically simulated results are in good agreement both with measured data and theoretical predictions.
While the experimental radargram (Figure 9a) was computed from real radar pulses distorted by the ionospheric phase dispersion, and therefore was corrected for it, in the simulated radargram shown on Figure 9b ionospheric phase shifts were not taken into account. To demonstrate the ionospheric correction algorithm performance, the results of simulations with accounting for ionospheric frequency dispersion of the phase and its partial correction is presented in Figure 9d separately. It can be seen that the simulated results match the measurements with a reasonable degree of agreement.
Some inaccuracies may arise from intrinsic errors of the Kirchhoff approximation applied here for the electromagnetic wave scattering description. On the other hand, they may occur due to some object features (variable material properties, etc.), completely missing in our simple model of perfectly reflecting surface of a given shape.
Ionospheric distortions, being properly compensated, also do not significantly impact the final result ( Figure 9e). To show how much distortion was in fact caused by the ionosphere in MEX orbit 9466, simulated radargram uncorrected for the ionospheric distortions is also shown in Figure 9c.
The ionospheric parameters chosen for the simulations correspond to the true conditions of the Martian atmosphere, which took place during the MEX 9466 orbit measurements. This orbit crosses the solar terminator in its northern part and in its southern part, which is analyzed here; the solar zenith angle changed within a range from 61 to 66 ∘ . According to radio occultation measurements (Zhang et al., 1990) and MARSIS results, the maximal plasma frequency of the Martian ionosphere exceeds 3.5 MHz for the solar zenith angles (SZA) smaller than 50 ∘ , approaching 4 MHz at SZA < 40 ∘ and exceeding this value during solar events. The NeMars ionospheric model predicts the maximal plasma frequency values about f p ≈ 3.4 MHz and f p = 4.14 MHz for the SZA = 60 ∘ and SZA = 0 ∘ , respectively. Thus, to make a complete test of the ionospheric correction algorithm performance, we investigate a hypothetical worst case scenario with SZA = 0 ∘ (the Sun is all the time at the local zenith). Simulated radargram, corresponding to this scenario, is shown in the Figure 9e. One can see notable degradation of the radargram quality; however, it is still readable and interpretable.
Effectiveness of the aperture synthesis largely agrees with theoretical estimates. It can be seen in Figure 10 that the aperture synthesis is able to suppress side clutter by 15-25 dB. Strong echo at t = 100 μs delay probably comes from transversal direction, which is not suppressed by the aperture synthesis.
The problem of correction of systematic phase distortion and random scattering of the radar signal acting simultaneously, has once been investigated for the diffraction on the ionospheric irregularities by Ilyushin (2008). From the algorithmic point of view, the problem effectively reduces to the finding of the correction phasẽ( ), corresponding to the maximal value of the function being optimized (i.e., intensity contrast). Under plausible assumptions about the ionospheric irregularities structure, it has been shown that when random scattering is present, the peak of the contrast function C 2 I is slightly flattened and widened but still remains unique and unambiguous. So, as a conclusion, random scattering does not significantly complicate the elimination of systematic regular phase distortion in the radar signal. However, no realistic model of kilometer-scale irregularities in the Martian ionosphere is known. Due to that, simplified assumptions have been made in that study.
The results of the present investigation largely agree with that previous study. As one can see from these two figures (11 and 12), both graphs demonstrate a clear unique maximum of the contrast function. In the case of the realistic surface roughness model, the maximum is somewhat lower and wider than that for ideally flat surface; however, it can still be easily found algorithmically with standard optimization routines. Thus, the phase correction algorithm based on the contrast optimization principle proves to be stable and robust with respect to the surface clutter, as well as to the ionospheric random scattering.

Conclusions and Remarks
In the paper presented here, a clutter simulating program CLUSIM, capable of emulating realistic radar echoes coming from rough planetary surface, is developed and described in detail. A model of surface landscape topography, based on interpolation of Mars Orbiter Laser Altimeter data set, has been elaborated. Surface echo simulator, exploiting Kirchhoff approximation and capable of simulation of aperture synthesis, developed as a C++ programming language procedure suitable for parallel computing systems using OpenMP parallel programming standard. A sample synthetic radargram is computed and compared to realistic radar instrument data obtained by the MARSIS experiment on board Mars Express interplanetary space mission. Previously elaborated routine for ionospheric distortion correction is also tested with simulated radar echo together with vertical plasma density profiles provided by NeMars ionospheric model, built with data recorded with the other model of the MARSIS radar. Comparative numerical tests confirmed the robustness of the ionospheric phase distortion correction algorithm with respect to the surface clutter. From the computational point of view, the newest point of the presented study is the aperture synthesis simulation technique. The novelty in the ionospheric phase correction techniques is the assessment of side clutter impact on the correction algorithm performance.