Radar Imaging with EISCAT 3D

A new incoherent scatter radar called EISCAT 3D is being constructed in Northern Scandinavia. It will have the capability of producing volumetric images of ionospheric plasma parameters using aperture synthesis radar imaging. This study uses the current design of EISCAT 3D to explore the theoretical radar imaging performance when imaging electron density in the E region and compares numerical techniques that could be used in practice. Of all imaging algorithms surveyed, the singular value decomposition with regularization gave the best results and was also found to be the most computationally efficient. The 5 estimated imaging performance indicates that the radar will be capable of detecting features down to approximately 90x90 m at a height of 100 km, which corresponds to a ∼ 0.05◦ angular resolution. The temporal resolution is dependent on the signal-to-noise ratio and range resolution. The signal-to-noise ratio calculations indicate that high resolution imaging of auroral precipitation is feasible. For example, with a range resolution of 1500 m, a time resolution of 10 seconds, and an electron density of 2 · 1011m−3, the correlation function estimates for radar scatter from the E-region can be measured with an uncertainty of 5 10 %. At a time resolution of 10 s and an image resolution of 90x90 m, the relative estimation error standard deviation of the image intensity is 10 %. Dividing the transmitting array into multiple independent transmitters to get at multiple-input-multiple-output (MIMO) interferometer system is also studied and this technique is found to increase imaging performance through improved visibility coverage. Although this reduces the signal-to-noise ratio, MIMO has successfully been applied to image strong radar echoes as meteors and polar mesospheric summer echoes. Use of the MIMO technique for ISR should be investigated further. 15 Copyright statement.


Introduction
One of the measurement challenges in the study of the Earth's ionized upper atmosphere using incoherent scatter radars (ISR) is that the measurements often do not match the intrinsic horizontal resolution of the physical phenomena that is being studied.
Conventional ISR measurements are ultimately limited in the transverse beam axis direction by the beam width of the radar 20 antenna, which is determined by the diffraction pattern of the antenna. Even for large antennas, the beam width is typically around 1 • . The mismatch between geophysical feature scales and horizontal resolution obtained by a typical ISR antenna is demonstrated in Figure 1, which shows an image of auroral airglow taken in the magnetic field aligned direction. Overlayed on the image are the antenna beam diameters of three incoherent scatter radar antennas: EISCAT UHF, EISCAT 3D, and Arecibo.
It is clear that the auroral precipitation has appreciable structure on scales smaller than the beam size. A conventional ISR measurement in this case will provide plasma parameters that are averaged over the area of the radar beam, preventing the observation of small sub-beamwidth scale structure. Only a radar with an antenna size of the Arecibo Observatory dish (305 m) would provide an antenna beam width that approaches the scale size of auroral precipitation.
Another measurement challenge for ISRs is temporal sampling of the spatial region of interest. Single dish radar systems 5 can only measure in one direction at any given time, and the ability to move the beam into another direction depends on the speed at which the antenna can be steered. In addition comes the minimum integration time required to measure one position.
It takes a long time to sample a large horizontal region, and even then, the measurements of different horizontal positions are obtained at different times.
In order to increase the spatial resolution of radio measurements without resorting to constructing an extremely large con- 10 tinuous antenna structure, a technique called aperture synthesis imaging can be used (e.g., Junklewitz et al., 2016). It relies on a sparse array of antennas to estimate a radio image with horizontal resolution equivalent to that of a large antenna. The correlation between the received signals can be used to produce an image of the brightness distribution of the radio source.
This technique is widely used in radio astronomy to image the intensity of radio waves originating from different sky positions.
The application of aperture synthesis imaging for radar, i.e., aperture synthesis radar imaging (ASRI), has been used in 15 space physics for observing high signal to noise ratio targets (Hysell et al., 2009;Chau et al., 2019). There is a good amount of literature on ASRI techniques in two dimensions (range and one transverse beam axis direction) for imaging field aligned irregularities (e.g. Hysell and Chau, 2012, and references therein). There have also been several researches on imaging of atmospheric and ionospheric features in three dimensions, e.g. Urco et al. (2019) who applied it to observations of PMSE with the Middle atmosphere ALOMAR radar system (MAARSY), Palmer et al. (1998) on the middle and upper atmosphere 20 radar in Japan, Yu et al. (2000) with a simulation study, and Chau and Woodman (2001) on the atmosphere over Jicamarca.
The currently available horizontal resolution of ASRI is around 0.5 • with Jicamarca, but down to 0.1 • for strong backscatter (Hysell and Chau, 2012) in the case of field-aligned ionospheric irregularities; and 0.6 • with MAARSY for polar mesospheric summer echoes (PMSE) .
However directly on incoherent scatter in three dimensions there is little literature, but some approaches have been made, 25 like Schlatter et al. (2015), who used the EISCAT Aperture Synthesis Imaging array and the EISCAT Svalbard radar to image the horizontal structure of Naturally Enhanced Ion Acoustic Lines (NEIALs) and Semeter et al. (2009) who interpolated sparse independent PFISR-measurements to estimate the electron density variation over a 65 × 60 km area during an auroral event.
In radar imaging, the measurements are in the so called visibility domain. Ensemble averages of the cross-correlation of complex voltages between two antennas represents a single sample of the visibility (Woodman, 1997;Urco et al., 2018).

30
Throughout this article we will use "far field" for the region further away than the Fraunhofer limit of the radar. The region closer than the Fraunhofer limit we will call "near field". If the radar target is in the far field of the radar, the visibility domain is related via a Fourier transform to the horizontal brightness distribution or the radio image. Then, the measurements are samples of the Fourier transform of the spatial variation of backscatter strength, or brightness distribution, of the target (Woodman, 1997). The measurements are used to calculate the brightness distribution, or image, of the target. So far, most of the incoherent scatter radar imaging has been done with a single transmitter and multiple receivers, thereby using a single-input multiple-output (SIMO) system. The number of measurements and degrees of freedom is here determined by the number of receivers and their relative locations. Instead of using only one transmitter, multiple transmitters can be used when performing radar imaging. This allows increasing the number of visibilities that can be measured, which can result in improved imaging performance, as long as the signal-to-noise ratio is sufficiently high. This technique is called multiple-5 input multiple-output (MIMO) radar (Fishler et al., 2006). The MIMO technique for increasing spatial resolution has recently been demonstrated with the Jicamarca radar when imaging equatorial electrojet echoes (Urco et al., 2018) and also with the MAARSY radar for imaging PMSE . The primary technical challenge with MIMO radar separating scattering corresponding to multiple transmitters on receive.

EISCAT 3D EISCAT UHF Arecibo
EISCAT 3D, from hereon referred to as E3D, is a new multi-static incoherent scatter radar that is being built in Norway, 10 Sweden, and Finland (McCrea et al., 2015;Kero et al., 2019). The core transmit and receive antenna array will be located in Skibotn,Norway (69.340 • N,20.313 • E). There will be additional bi-static receiver antenna sites in Kaiseniemi,Sweden (68.267 • N,19.448 • E) and Karesuvanto,Finland (68.463 • N,22.458 • E). The core array of E3D will consist of 109 subarrays, each containing 91 antennas. The one-way half power full beamwidth (HPBW) or illuminated angle of the core array will be 1 • . On transmission, the array is capable of transmitting up to 5 MW of peak power at a frequency of 233 MHz. Additionally, 15 there are 10 receive-only outrigger antennas around the core array, providing longer antenna spacings that can be used for high resolution ASRI. Imaging will already be necessary to maintain the perpendicular resolution constant in the transition from EISCAT VHF and UHF to E3D. It is possible that the EISCAT 3D radar can also be configured as a MIMO system, where the core array is separated into smaller subarrays, which act as independent transmitters at slightly different locations. During the design phase, Lehtinen (2014) investigated the imaging performance of possible layouts of E3D in the far-field. The study 20 however does not include the current layout that is being built. EISCAT 3D will not be able to measure radar echoes from magnetic field aligned irregularities, so it will not be possible to assume that the scattering originates from a two dimensional plane where the radar scattering wave vector is perpendicular to the magnetic field. All radar imaging will need to be done in 3D and mostly for incoherent scatter. This poses several two main challenges: 1) the signal-to-noise ratio will in typical cases be determined by incoherent scatter, which is much smaller than 25 that used conventionally for ASRI; 2) there are more unknowns that need to be estimated, as at each range there is a 2D image instead of a 1D image that needs to be estimated. For E3D, the Fraunhofer limit is at 2D 2 /λ ≈ 2000 km, where D ≈ 1.2 km is the longest baseline and λ = 1.3 m is the wavelength of the radar. Measurements of the ionosphere are therefore taken in the nearfield of the radar. Woodman (1997) describes a technique to correct for the curvature in the backscattered field with an analogy of lens focusing. In this study, a 30 different approach has been taken, where the nearfield geometry is directly included in the forward model of the linear inverse problems formalism. In this case, it is not possible to resort to frequency domain methods to diagonalize the forward model.
This comes at an increase in computational complexity, but is not prohibitive in terms of computational cost with modern computers.
In this study, we will simulate radar imaging measurement capabilities of the upcoming EISCAT 3D radar. The study is divided into the following sections. In Sect. 2, we investigate the achievable time and range resolution of E3D, and how they are connected. An expression for the cross-correlations between the received signals, taking into account the nearfield geometry, is derived in Sect. 3. Sect. 4 describes the nearfield forward model for radar imaging and describes several numerical techniques for solving the linear inverse problem. This section also includes a study of imaging resolution based on simulated 5 imaging measurements.

Time resolution
In this section, we will calculate the required integration time for a certain range resolution with E3D. The elementary radar imaging measurement is an estimate of the cross-correlation of the scattered complex voltage measured by two antenna modules. The integration time in this case is the minimum amount of time that is needed to obtain a error standard deviation for the 10 cross-correlation estimate that is equal to a predefined limit. The estimation error of the cross-correlation determines the measurement error for the imaging inverse problem. By investigating the variance of the cross-correlation estimate using statistical properties of the incoherent scatter signal, we can decouple the problem of time and range resolution from imaging resolution, allowing us to study the performance of the imaging algorithm with a certain measurement error standard deviation.
Our signal-to-noise calculations will be based on an observation of incoherent scatter from ionospheric plasma, which is 15 the case with the smallest expected signal-to-noise ratio. We have ignored self-clutter, as the combination of the E3D core transmitter illuminating the target and a single receiving subarray receiver module will inevitably be within a low signal-tonoise ratio regime, dominated by receiver noise.
We will first deduce an expression for the measurement rate, that is, how many measurements are taken per second. There are two factors that determine the maximum rate at which independent observations of the scattering from the ionosphere can 20 be made: 1) The minimum inter-pulse period length, which we set to d/τ p , with d the duty-cycle and τ p the pulse length. 2) The incoherent scatter decorrelation time, which is inversely proportional to the bandwidth of the incoherent scatter radar spectrum B. The maximum of these two time scales determines the frequency of independent measurements that can be made: If a transmitted longpulse is divided into N P bits, the number of measurements per longpulse can be multiplied by N P . In 25 the E region, we can assume that the autocorrelation function is constant for the purpose of estimating the variance. Then the number of lagged product measurements per transmit pulse is N P (N P − 1)/2 because we also can use measurements with different time lags. For sake of simplicity, we assume that all lags within a radar transmit pulse are equally informative. This is approximately the case for E-region plasma measured using E3D. The number of measurements per second is then

30
Next, we will estimate the number of measurements needed for reducing the measurement error of an average crosscorrelation measurement to a certain level. We consider a measurement model where a measurement m is described by a linear combination of the parameter we want to estimate m = x + ξ, where x and ξ are considered as proper complex Gaussian random variables with zero mean and variance of respectively P S and P N . The noise power estimate P N is assumed to have no error. We estimate the signal power witĥ where the bar denotes complex conjugation. It can then be shown that where ε is the relative standard deviation and K is the number of measurements (Farley, 1969). If we require the correlation function to have a relative uncertainty under a certain level , f.ex. = 0.05 = 5%, the equation can be solved for K in order to get the number of needed samples

10
where SNR is the signal-to-noise ratio SNR. The integration time required to obtain a measurement with a certain level of uncertainty is now or written as a function of SNR 15 The received signal power P S can be found by the radar equation where P tx is the transmitted power, G tx is the transmit and G rx the receive gain, λ is the radar wavelength, σ is the scattering cross-section and R tx and R rx are the distance from the scattering volume to the transmitter and receiver (e.g. Sato, 1989).
Assuming that the Debye length is much smaller than the radar wavelength, the effective scattering cross-section for a single 20 electron in plasma (Beynon and Williams, 1978) is Here, σ e is the Thomson scattering cross-section σ e = 4πr 2 e , T i is the ion and T e the electron temperature. The total scattering cross-section can be found by adding up the cross-sections of all electrons in the illuminated volume The scattering volume can be approximated using a spherical cone: Here, ∆r = cτ b /2 is the range resolution of the measurement, where τ b = τ p /N p is the baud length, r is the range of the volume and θ is the HPBW angle of the radar. By range, we mean the range from the center of the core array in Skibotn to the target. 5 We assume that the noise is constant through the ion line spectrum. The noise power is then given by where T sys is the system noise temperature, B is the bandwidth of the incoming ion-line, and k B is the Boltzmann constant.
The bandwidth we assume to be equal two times the ion thermal velocity times wave number (2v th k) or the inverse of the pulse length τ −1 b , depending on which one is the largest. The ion thermal velocity is given by where m i is the ion mass, which we set equal to 31 u corresponding to a mixture of O + 2 and N O + . These are the two most dominant ion species in the E region (Brekke, 2013). The system noise temperature we set to 100 K. For all bit lengths investigated here, τ −1 b exceeds 2v th k by at least one order of magnitude. The bandwidth is therefore independent on the ion composition as long as the measurements are restricted to the E region. 15 We can now use this to calculate the integration time of electron density measurements in the E-layer at 150 km for different range resolutions. The standard deviation requirement is set at 5%. The electron and ion temperatures we set to 400 K and 300 K, respectively. We assume a monostatic radar with frequency f = 230 MHz, HPBW θ =1°, and a transmitter with power of 5 MW. The transmitter gain for the core array we set to 43 dB and the receiver gain to 22 dB for one subarray of the imaging array. The interpulse period τ IPP is 2 ms and the longpulse length is 0.5 ms. The results are shown in Fig. 2. 20 The figure shows that the integration time decreases with increasing electron density and decreasing range resolution. This confirms the expected tradeoff between range and time resolution. If the electron density is not too low, a time resolution of a few seconds is possible. This however assumes a relatively low range resolution of 1000 -2000 m, which still provides some useful information about the E-region plasma. When keeping a constant standard deviation, an enhanced electron density can be used to improve either the time or range resolution.

25
When using MIMO imaging, the core array is divided into multiple independent groups when transmitting. This provides more baselines as well as increases the maximum antenna separation. In this case, the imaging resolution will be improved by having a larger aperture. One of the challenges in this case will be separating the signals from different transmitters in the case of overspread radar targets. We assume that separate transmitters operate at the same frequency and that the transmitted signals are distinguished using radar transmit coding. This can be achieved in practice using a different pseudorandom transmit code smaller antenna area, the transmit gain must be divided by the number of transmitters. Additionally, there could be crosscoupling between antennas, which might cause buffer zones between transmitters. Then the antenna area and gain decrease furthermore. In conclusion, the integration time for MIMO will at least be the number of transmitters times the integration time for SIMO.
The calculations do not include other echoes than from incoherent scatter or other enhancements than electron density. In the 5 case of PMSE (e.g. Urco et al., 2019) and NEIALS (e.g. Grydeland et al., 2004;Schlatter et al., 2015), the echo is significantly stronger than for incoherent scatter. These enhancements will also make shorter integration time available, and will be more promising candidates for use of MIMO imaging.

Baseline cross-correlation
In this section we calculate the correlation between signals from two different baselines, that are transmitter-receiver-pairs. The aim is to determine which baselines provide information about ionospheric features of a certain scale size, and to determine how the nearfield geometry affects this correlation.
We consider a case with one transmitter and two receivers placed with a equally long distance from the transmitter in every 5 direction. This configuration is shown in Fig. 3. Let the transmitter be placed in the origin and the receivers at |P 1 and |P 2 .
Let the transmitter transmit an signal of the form V 0 = Ke iωt , where K is a time-independent constant, ω, is the transmit frequency and t is time. The electrical potential induced to receiver antenna r then becomes where T in = ||R i +r n |/c is the time delay from the transmitter to scatterer n and T srn = ||R sr −r n |/c is the time delay from 10 scatterer n to receiver r. Here |R i is the vector from the transmitter to the centre of the illuminated plasma volume, |R sr is the vector from the centre of the plasma volume to receiver r, and |r n is the vector from the centre of the plasma volume to scatterer n, like in Fig. 3. N is the number of scatterers in the scattering volume, and G ∈ is the scattering gain which includes the free-space path loss. The gain may be dependent on the position of the scatterer, the scatterer itself, and on time, but we neglect these dependencies. We also neglect that the distance to the scatterer varies between the transmitter-receiver baselines. This has a order of magnitude of~10 m, which is lower than the best available range resolution.
The cross-correlation function for time lag τ = 0 can then be written as By taking the first-order Taylor approximation of the time delays around |r n = |0 , we get that where the hat denotes a unit vector. Carrying out this approximation is essentially the same as assuming plane waves. When keeping also the second order terms, the nearfield correction described by Woodman (1997) can be deduced. 10 We note that − ω We assume that the scatterer positions are independent identical normally distributed with mean |µ and covariance L, that is like a Gaussian blob, 15 We use the definition of expectation and then solve the integral. Since the positions of the scatterers are assumed to be independent, the expectation becomes zero when n = n . In addition, in first order approximation, R s2 − R s1 ≈ 0 because D << h.

The result then becomes
The normalized cross-correlation function, We note that if the transmitter(s) and all receivers lay in a plane, the vertical components of the Bragg scattering vectors are exactly equal and make the vertical components of |µ and L, namely µ z , L xz , L yz and L zz arbitrary. This means that the 25 horizontal resolution is independent on the vertical resolution.
Equation 20 for µ| = [0, 0, 0] and L = (L/2) 2 I is plotted in Fig. 4, where L is the extent of the ionospheric feature in all dimensions, and I is the identity matrix. Figure 4 also shows numerically simulated normalized correlation based on a direct simulation of Equation 16, which does not significantly differ from the analytical expression. The plot shows that for a height of 10 5 m (100 km) and a baseline of 211 m, the correlation crosses 0.95 at a blob size of 70 m and 0.5 at a blob size of 250 m. At 100 km height, the radar beam of E3D is about 1800 m wide. This means that when considering a maximum baseline 5 of 200 m and a ionospheric feature that is larger than 250x250 m, addition of longer baselines contribute less to recover the image. The E3D core has a maximum baseline of 75 m. We can simplify these calculations by setting the magnitude of the desired least cross-correlation to R, We assume that the scatterers have equal variance in x-and y-direction (L xx = L yy = (L/2) 2 ) and that all directions are 10 uncorrelated (L xy = L xz = L yz = 0). By using the geometry as in Fig. 3, By combining Eqs. 20, 21, and 22, we get , which can be rewritten to an expression for the feature size L: For short baselines or long distances D h/5, the expression can be simplified, we solve for the baseline D and get The resulting expression shows how long the baseline can be to still get contribution to recovering the feature. Equation 24 is plotted in Fig. 5.

20
A longer baseline can contribute to recover smaller features, but the improvement will decrease the longer the baseline is.
For example if we want to resolve a feature with size 100 m, baselines up to 200 m have large contribution to the imaging.
Adding longer baselines will improve the resolving less, and stopping slightly above 1 km . This means that the improvement of the imaging quality by including the E3D outriggers will be large for the closest outriggers. The signal received by furthest ones will correlate little with the signal received by the core. From Fig. 5 we see that the correlation in the longest E3D baseline  When inserting R = 0.01, Eq. (24) shows that the diffraction limit is the same as for planar scatter under the assumptions mentioned in the deduction.
Baselines between the receiver sites in Skibotn, Karesuvanto and Kaiseniemi are so long that they can not be used for imaging as signals will not be correlated anymore. The baseline cross-correlation calculations also do not claim that the image is well recovered if including the largest baseline. This is more dependent on which baselines are used, how they are distributed,

Radar imaging model
We consider a radar that may have single or multiple inputs (transmitters), and multiple outputs (receivers) (SIMO or MIMO).
The radar illuminates a plasma volume at range R with thickness dr, and inside of the one-way HPBW θ. We imagine that the volume is divided into M parts, or pixels, see Fig. 6 .
The signal transmitted from transmitter A and spread by plasma element/pixel q causes a voltage fluctuation in the receivers. 5 The voltage fluctuation of receiver D due to transmitter A and plasma pixel q is denoted as V q AD = F V A e 2πif T q AD , where V A is the amplitude of the signal sent by transmitter A, F is a function of the received signal amplitude, f is the radar transmitting frequency and T q AD is the time delay of the signal due to travelling from transmitter A via pixel q to receiver D, c.f. Fig. 6. The correlation between signals from two different baselines AD and BE due to an infinitesimal scattering volume dV can be described as where P t is transmit power, G t is transmit gain, G r is receiver gain, λ is the radar wavelength, σ p is the scattering cross-section for a single electron given by Eq. (9), n e is the electron density, R t is the distance from transmitter to the scattering volume, R r is the distance from the scattering volume to the receiver, and |r is the position of the scattering volume. We integrate over the 5 whole scattering volume to get the whole measurement. At a certain time lag, we get the correlation for the range of interest.
We assume that the gains are constant inside of the radar beam and zero otherwise and neglect the dependency of R r and R t on the exact position of the scattering volume. The correlation can then be written as We assume that the electron density (or brightness) distribution can be written as a sum of discretized parts with constant electron density. We neglect variations in the phase shift inside of one part. The integral can then be replaced with a sum The first factor here is constant and can be normalized away. The number of discretizations Q is still needed in the simulations if the original image has an other resolution than the reconstructions. The series of measurements can be written on matrix 5 form, is the theory matrix, and |ε = [ε AAAA , ε AAAB , . . . , ε KKKK ] H is the complex normally distributed noise vector. 10 Sometimes it is more convenient to have the crosscorrelations on matrix form. The measurements can be transferred from the one form to the other simply through reshaping the vector |m to a matrix M or opposite: To get an estimate of the intensities of the plasma in the image, Eq. (28) has to be inverted so that 15 where B is a matrix that reconstructs the image |x from the measurements |m . When inserting Eq. (28) into Eq. (30) and neglecting noise, we get |x = BA|x . We wish that the reconstructed image is as close to the reality as possible, and so taking B = A −1 would give a perfect solution. However, since we have an underdetermined problem, A cannot be directly invertible.
Other attempts are therefore needed.

Matched filter 20
When the scatterers are behind the Fraunhofer limit in the far field, Eq. (28) represents a Fourier transform. One approach to get back the original image would be the inverse Fourier transform, which can be represented the hermitian conjugate of the theory matrix, B = A H like a matched filter (MF). Unfortunately, the samples of the Fourier-transformed image, that are the visibilities, are sparsely and incomplete scattered and the problem gets underdetermined (Hysell and Chau, 2012;Harding and Milla, 2013). The approach can be interpreted as steering the beam after the statistical averaging and is therefore also called beamforming.

Capon method
Another approach is the Capon method (Palmer et al., 1998). The purpose of this method is to minimize the intensities in all 5 other directions than the direction of interest, that is to minimize the sidelobes of the antenna array in directions with interfering sources. The result is to invert the matrix of correlation measurements M (Palmer et al., 1998). In order to continue using the notation in this article, M −1 is reshaped back to a vector |m −1 . The estimated intensities by the Capon method can then be written as 10 where the fraction denotes element-wise division.

Singular value decomposition
The problem in Eq. (28) is overdetermined if the number of unknowns, that is the number of discretizations, is less than the number of measurements. This can be the case if we solve for a imaging resolution that is low enough. We can then use the method of least squares to solve it, getting One can also use the singular value decomposition (SVD) on the theory matrix, A = USV H , where S is a diagonal matrix containing the singular values, that are square roots of the eigenvalues of A H A, V contains the normalized eigenvectors of A H A, and U contains the normalized eigenvectors of AA H . The inversion matrix B can then be written as 20 which can be shown still gives the same solution as ordinary least squares, but with increased numerical accuracy (Aster et al., 2013). Because of inverting the singular values, the eigenvectors corresponding to the smallest values contribute most to the variance of the solution and make the solution sensitive to noise. Also, the problem can be rank deficient, that is that several columns in the theory matrix are nearly linear dependent on each other. The problem is then said to be ill-conditioned or multicollinear.

25
In such cases, some singular values will be practically zero and the solution may be hidden in the noise. To prevent the noise sensitivity, the solutions can be regularized. This makes the reconstruction biased towards smoothness and zero, but less noisy (Aster et al., 2013). We here consider two regularization techniques, truncated SVD (TSVD) and Tikhonov regularization. In TSVD, the inverse of the singular values below some limit are set to zero. The eigenvectors corresponding to the smallest singular values will then not contribute to the result. These eigenvectors often contain high frequency components. Ignoring 30 them makes the solution smoother. Tikhonov regularization or ridge regression can be done in several ways. In this article, we use zeroth order Tikhonov regularization where the singular values s i are inverted with where α is a regularization parameter. By using SVD, we also can get the variance |Σx of the estimates. For pure least squares, it is diag((A H A) −1 ), for regularized least squares it is

CLEAN
CLEAN is another attempt to reduce sidelobes. It is based on the matched filter approach, but iteratively finds the real structure in the field of view (Högbom, 1974). It supposes a source where the image reconstructed by the matched filter is brightest.
The source is added to an image that only is containing the suspected sources, which will be the reconstructed image. Then 10 the measurements that the radar would have measured if the reconstructed image were the true image are subtracted from the real measurements and the next suspected source is found. This procedure is repeated until there are no clear sources left in the measurements (Högbom, 1974). The method is a special case of compressed sensing requires an assumption on how the measured sources look like (Harding and Milla, 2013). For sparse sources, a Dirac delta function could be appropriate, but lead to sparse solutions.

Performance of the radar layouts
We considered different radar layouts. The layouts together with plots of the visibilities and the point spread function are shown in Fig. 7.
When considering a layout with multiple transmitters and multiple receivers (MIMO), it is assumed that the signals from different transmitters can be distinguished. This increases the number of virtual receivers and thereby the visibilities get more 20 widespread and denser, cf. Fig 7. However, using multiple transmitters increases the integration time as described in Sect. 2.
We note that when receiving with the outriggers, the main beam becomes narrower. Also, there are gaps in the visibilities. This is due to the sparse locations of the outriggers and makes the point spread function look more irregularly (cmp . Fig 7c and 7f).
With multiple transmitters, the main beam becomes even narrower (cmp . Fig 7c and 7i). When both using multiple transmitters and receiving outriggers, the gaps in the visibility domain partially get filled and the sidelobes are clearly reduced. The MIMO 25 layout used here could possibly be improved by using positions of the transmitters so that gaps in the visibility get more filled.

Performance of the imaging techniques
We simulate E3D measurements using Eq. 28 and with the presented antenna configurations. As original image, we use a part of Fig. 1 were truncated. This value gives a good compromise of resolution and low noise level (regularization). For CLEAN, we used a gain of 1 and a threshold of 1.36 times the average value. We tried both Dirac delta and Gaussian functions in the CLEAN kernel. In capon filtering, it happens that the correlation measurement matrix M is singular. In such cases, TSVD is used to invert the matrix. This truncation ignores singular values that are less than 0.03 % of the largest singular value. 5 Noise is added to all cross-correlations which corresponds to white complex Gaussian noise with zero mean and 5 % standard deviation in each receiver. The noise is equal for every reconstruction of a single resolution, but varies between reconstruction in different resolutions. The results for the SIMO layout is shown in Fig. 8.
Of the reconstruction techniques, TSVD clearly gives the best results. It is also the only method that fairly reproduces the shape of the true image. Capon also partly reproduces the shape, but far worse than TSVD. The matched filter apparently only 10 reproduces something similar to the point spread function. The performance of CLEAN (not shown here) is accordingly poor.
In terms of calculation time, CLEAN is the slowest algorithm followed by TSVD. MF and Capon are relatively fast. These differences get stronger when also considering MIMO. Most of the computation time of TSVD is used to invert the theory matrix. Since the theory matrix only varies from experiment to experiment, it must only be inverted once and can be saved afterwards. The computation time therefore is reduced to a simple matrix multiplication, and it is not considered as a problem 15 for the real radar. We therefore concentrate on images reconstructed with techniques using SVD. Here, we compared ordinary least squares, TSVD with truncating singular values under 2 %, like before, and Tikhonov regularization with regularizing parameter α = 10 and 100. These results are shown in Figs. 9 to 12 for the four layouts shown in Fig. 7, that are SIMO without and with outriggers and one MIMO case without and with outriggers, respectively.
For all layouts, at a considered resolution of 20x20 pixels in the radar main beam (here, 1 pixel≈100x100 m), the image 20 reconstruction with method of least squares is very noisy. The singular values of A are varying over several orders of magnitude which is a sign of that columns in A are linear dependent on each other. The regularized solutions look considerably better, with a regularization parameter of 10, the recovered images are still a bit noisy, but with stronger regularization the images get smoother and closer to the original.
The two most strongly regularized images to the upper-right and bottom-middle contain stripes if the radar layout is including 25 the outriggers. This is propably because the visibility in some regions has gaps, cf. Fig 7. When only considering the core array, there are no gaps other than the spacing between antennas. The recovered images without the outriggers look smoother than including the outriggers, but when including the outriggers, more details of the original image can be seen. Also, in the MIMO case with outriggers, the feature in the south-east can be seen in the reconstruction. For the other layouts it is less visible and not clearly distinguishable from the main feature in the north.

30
The uncertainty of the reconstruction itself is given by the variance of the recovered image, Eq (34). The mean standard deviation for the different layouts and reconstruction techniques is shown in Fig. 13. The plots of the least square variance are comparable to the variance plots in Lehtinen (2014). We note that while Lehtinen (2014) investigates far-field imaging, Fig. 13 shows near-field imaging.  By using the standard deviation we neglect errors introduced by the discretization because they are not included in the variance. This assumption is true if the true image has the same resolution as the reconstruction, but that is only for the case of what Kaipio and Somersalo (2010) call a "inverse crime". In reality, the target of E3D, the electron fluctuations in the ionosphere, is not discrete with steps of several metres. Also, by regularization, bias is introduced to the solution, which the variance does not take into account. Therefore, we also used the similarity to the true image for uncertainty estimation. As 5 a measure, we used the mean square deviation, s = N i=1 (xi−xi) 2 N so a low value of s means great similarity. Because the original image and the reconstruction have different resolutions, the smallest is scaled up. The scaling was done by Lanzcos resampling with a cos 2 -kernel. A drawback with the MSE is that it could be influenced by the target, while the variance is not.
The mean standard deviation and the similarity to the original image are shown in Fig. 13 for all layouts considered here and reconstruction resolutions up to 100x100 pixels. Figure 10. Comparison of reconstructions using SVD for the SIMO case including the outriggers. Else, the plots corresponds to Fig. 9. Figure 11. Comparison of reconstructions using SVD for the MIMO case only using the core antennas. Else, the plots corresponds to Fig. 9. The variance of the recovered image is strongly increasing with the resolution we assume/want that the radar target has.
Image recovering with LS gives the highest variance for all layouts. The reason is the multicollinearity in the theory matrix which amplifies the noise in the recovered image. At small resolution, the variances are equal for the different reconstruction techniques, but diverge when the regularization starts to influence the results. This divergence happens later when including the outrigger antennas and also later when using MIMO than for SIMO. For high resolutions, the variance of the TSVD solution is 5 the lowest. However, since bias is introduced by regularizing the solution, this does not necessarily mean that TSVD solution is the best.
The mean square error (MSE) of the recovered images is in general higher than their mean variance. For small resolutions, it decreases with increasing resolution until it reaches a bottom point. The error then increases again. For LS and Tikhonov with α = 10, the minimum is at 10-20 pixels per direction. When including the outriggers, the minimum is at a later stage. Also, 10 the error is lower. We also note the dip of error at 97x97 pixels. This is exactly where the resolution of the recovered image matches the resolution of the original image so these dips are the effect of inverse crimes and therefore not transferable to the real radar. For high resolutions, the MSE is higher for MIMO than for SIMO when using Tikhonov. This could indicate that for MIMO, more regularization is required.
The original image contains values between 0 and 1 m -3 with a mean of about 0.5 m -3 . In real, the values will be far higher 15 and the uncertainty will increase accordingly. Therefore, the standard deviation and the MSE are plotted relative to the mean value of the original image. In order to have a good recover, the relative mean error should be below 1 and, if possible, far below that. All regularized solutions would fulfil this criterion, but the two strongest regularizations clearly have the lowest MSE. The minimum of MSE seems to be somewhere between 60x60 pixels for MIMO and 90x90 pixels for SIMO. In practice, the image reconstructions for higher resolutions look very similar to low resolution (20x20) without adding more details but with better quality of the reconstructed image.
In the MSE plots, the curves flatten out to a minimum relative MSE at about 10 %. At 100 km range, 20x20 pixels corre-5 sponds to a resolution of around 90x90 m. The TSVD indicates that the recovered image with MIMO could be improved with stronger Tikhonov regularization, but this has not been investigated.
The MSE of TSVD does not decrease significantly from SIMO to MIMO. Therefore it seems that there is little gain by using the MIMO layouts considered here. However, the feature in the bottom right part of the image in Figs. 10 and 12 gets clearer with MIMO. For other targets, these results may look different. When comparing the point spread functions in Fig. 7, it could 10 be that MIMO configuration is better for point-like targets, like space debris or meteors, but this is beyond the scope of this article.
In this article, the MIMO approach to ISR and E3D has only been treated superficially. There are still some questions that must be answered. To distinguish between the transmitters, we assumed code diversity. However, there is need to study how well the signals can be distinguished. This can influence the possible number of transmitters. The placing of the transmitters 15 was not investigated here, the example in this article is a simple proposal. At the same time, the transmitter locations can have great influence on the visibility coverage. Also, the SNR and integration time calculations for MIMO would need to be investigated more thoroughly.

Conclusion
In this article, we have studied the temporal and spatial resolution of the upcoming E3D radar in the case of aperture synthesis 20 radar imaging, primarily focusing on the feasibility of imaging the incoherent scatter radar return from the E-region. The most up to date radar design specifications at the time of writing this article was used as a basis of this study.
We find that the range and time resolutions are dependent on each other. When keeping the uncertainty level constant, a better range resolution goes on the cost of the time resolution. With an increase in the electron density, the resolution in time and/or range can be improved without increasing the noise level. Under normal conditions in the E-layer (T e ≈ 400 K, T i ≈ 300 K, 25 n e ≈ 10 11 m -3 ), with a desired integration time of 10 seconds, the achievable range resolution is slightly more above 1500 m.
The horizontal (imaging) resolution depends on the radar layout and the imaging technique. The imaging techniques that were evaluated were: matched filter, least squares using singular value decomposition without and with regularization, Capon, and CLEAN. Of these techniques, only regularized least squares gave satisfactory results. The two regularization techniques of either truncating or damping of the inverse singular values both worked and gave similar results.

30
These image reconstructions can be reduced to a simple matrix multiplication by saving the inverted theory matrix. Regularized SVD is therefore is among the fastest reconstruction techniques amongst the ones evaluated. With Tikhonov regularization with a damping coefficient of 100, or truncating away singular values below 2 % of the largest value, the relative error of the recovered image can go down to 10 %. The resolution of the recovered image is about 60x60 pixels, at 100 km range this corresponds to 30x30 m, but features smaller than 90x90 m will be blurred out.
The simulation results show that using the outriggers increases the imaging accuracy. Dividing the core array into multiple transmitters to get a MIMO system seems to increase the imaging resolution if the target is smooth. MIMO also has the drawback that it needs stronger signals or more integration time to keep the same measurement accuracy as SIMO. However, 5 this needs further investigation. as MIMO may be useful for very bright targets such as PMSE, as well as point-like targets like space debris or meteors, but the latter needs further investigation.
We conclude with that radar imaging with EISCAT 3D is feasible.