Enhancing the spatiotemporal features of polar mesosphere summer echoes using coherent MIMO and radar imaging at MAARSY

Polar mesospheric summer echoes (PMSEs) are very strong radar echoes caused by the presence of ice particles, turbulence, and free electrons in the mesosphere over polar regions. For more than three decades, PMSEs have been used as natural tracers of the complicated atmospheric dynamics of this region. Neutral winds and turbulence parameters have been obtained assuming PMSE horizontal homogeneity on scales of tens of kilometers. Recent radar imaging studies have shown that PMSEs are not homogeneous on these scales and instead they are composed of kilometer-scale structures. In this paper, we present a technique that allows PMSE observations with unprecedented angular resolution (∼0.6^∘). The technique combines the concept of coherent MIMO (Multiple Input Multiple Output) and two high-resolution imaging techniques, i.e., Capon and maximum entropy (MaxEnt). The resulting resolution is evaluated by imaging specular meteor echoes. The gain in angular resolution compared to previous approaches using SIMO (Single Input Multiple Output) and Capon is at least a factor of 2; i.e., at 85 km, we obtain a horizontal resolution of ∼900 m. The advantage of the new technique is evaluated with two events of 3-D PMSE structures showing: (1) horizontal wavelengths of 8–10 km and periods of 4–7 min, drifting with the background wind, and (2) horizontal wavelengths of 12–16 km and periods of 15–20 min, not drifting with the background wind. Besides the advantages of the implemented technique, we discuss its current challenges, like the use of reduced power aperture and processing time, as well as the future opportunities for improving the understanding of the complex small-scale atmospheric dynamics behind PMSEs.

1 Introduction

The so-called MIMO (Multiple Input Multiple Output) technique is being widely used in the fields of telecommunications and radar remote sensing (e.g Telatar, 1999; Huang et al., 2011; Foschini and Gans, 1998). Recently Urco et al. (2018) have shown that the use of multiple transmitters and multiple receivers can significantly improve the angular resolution of coherent atmospheric and ionospheric radars. In that work, MIMO was used to observe equatorial electrojet (EEJ) field-aligned irregularities at Jicamarca in combination with the well-established radar imaging technique Capon (e.g., Palmer et al., 1998). The multiple transmitter part was implemented with three different diversity schemes, i.e., temporal, code, and polarization. The resulting angular resolution was superior, by at least a factor of 4, to previous efforts using a single transmitter and the same receiving configuration, i.e., SIMO (Single Input Multiple Output). Given that the EEJ irregularities are field-aligned with the Earth's magnetic field, angular imaging was performed only in the magnetic east–west direction.

Based on this successful implementation, we decided to implement coherent MIMO to improve the angular resolution of the Middle Atmosphere ALOMAR Radar System (MAARSY) (16.04^∘ E, 69.30^∘ N) and to study polar mesospheric summer echoes (PMSEs). PMSEs present strong radar cross sections (RCSs) that allow them to be observed with less transmitting power, which is the case when using MIMO. Previous efforts to study their spatial structure have been limited to a few kilometers' spatial resolution and a few minutes' temporal resolution (e.g., Yu et al., 2001; Latteck et al., 2012a; Stober et al., 2013; Sommer and Chau, 2016). Recently, Stober et al. (2018) have presented many examples of monochromatic gravity waves (GWs) and Kelvin–Helmholtz instabilities (KHIs) using 9 days of multibeam PMSE observations with MAARSY.

PMSEs are strong echoes, more than 50 dB stronger than expected echoes from free electrons in the D region, and there is a consensus that they are generated by atmospheric turbulence and require the presence of free electrons and charged ice particles (e.g., Rapp et al., 2002; Varney et al., 2011, and references therein). Although PMSEs have been studied since the late 1970s (e.g., Ecklund and Balsley, 1981; Hoppe et al., 1988; Kelley and Ulwick, 1988; Havnes et al., 1996; Rapp and Lübken, 2004), until recently they have been considered very aspect-sensitive and homogeneous on scales of a few tens of kilometers, at least when observed at very high frequencies (VHFs) (e.g., Czechowsky et al., 1988; Zecha et al., 2001; Yu et al., 2001).

Based on recent multibeam observations as well as radar imaging, Sommer and Chau (2016) have concluded that the PMSEs are not as aspect-sensitive as previously reported, and instead, most of the time they are organized in kilometer-scale spatial structures drifting across the observing beams. Such results have been independently verified with bistatic observations at VHFs, where PMSEs were observed with small systems at zenith angles close to 30^∘ (e.g., Chau et al., 2018).

The results of Sommer and Chau (2016) were obtained with MAARSY using the whole antenna array for transmitting and an antenna compression approach, i.e., a wide beam by properly phasing the antennas (e.g., Woodman and Chau, 2001) and a multiple-receiver configuration. The spatial structures were obtained using the Capon technique due to its implementation simplicity and its relatively fast processing speed.

Given that PMSEs are highly associated with noctilucent clouds (NLCs) (e.g., Hoppe et al., 1990; Stebel et al., 2000; Kaifler et al., 2011), spatial structures ranging from a few hundreds of meters to a few tens of kilometers observed in NLCs (e.g., Baumgarten and Fritts, 2014) are expected to be observed also in PMSEs. Indeed this is the case, PMSE structures of a few kilometers have been already reported by Sommer and Chau (2016) and structures of a few tens of kilometers have been reported by Chau et al. (2018).

Although progress has been made in discriminating between spatial and temporal ambiguities in PMSE observations, the achieved angular resolution has been mainly limited by two factors: (1) the effective area in the visibility plane and (2) the number of independent spatial samples (e.g., Woodman, 1997). By implementing MIMO, we are able to improve both, i.e., a larger effective area and a higher number of independent visibility samples. In addition, by implementing maximum entropy (MaxEnt), which is more computationally demanding than Capon, we are able to further improve the angular resolution (e.g., Hysell and Chau, 2006).

In this work, we have implemented coherent MIMO at MAARSY using 3 spatially separated antenna sections on transmission and 15 on reception. Moreover, time diversity was employed in order to isolate radar echoes corresponding to each transmitting section; i.e., the transmitters were interleaved every 4 ms. The resulting effective number of virtual receivers by using MIMO was 45 and the angular resolution achieved was ∼0.6^∘. It is equivalent to an antenna area of 450 m diameter, more than 5 times larger than the nominal diameter of the MAARSY antenna.

Our paper is organized as follows. We first present the experiment configuration with a specific emphasis on the MIMO implementation. Then we describe the radar imaging implementation for both Capon and MaxEnt techniques. The PMSE results are shown in Sect. 4 for SIMO and MIMO using both Capon and MaxEnt. Within this section, two events are studied in detail, one in which the observed waves drift with the background wind and a second one in which the waves do not propagate with the wind. Finally, the results of our MIMO implementation are discussed, followed by conclusions.

2 Experiment configuration

2.1 MAARSY

MAARSY is an active phased antenna array operating at 53 MHz, located in Andoya, Norway (69.30^∘ N, 16.04^∘ E). The array consists of 433 antenna elements, each with its own transceiver module that allow us to modulate the antennas in phase and amplitude independently. Using this capability, the transmitting or receiving beam can be steered in a desired direction up to 30^∘ off zenith, with an angular resolution of 3.6^∘ (e.g., Latteck et al., 2012a). In addition to its multibeam capability, MAARSY can be used for in-beam imaging experiments. In this case, the signals from a selected number of receiving antennas are stored, and later a digital beamforming algorithm (imaging) is applied to the data. Unlike the multibeam experiment, imaging allows a 2-D image to be obtained at once, avoiding the interleave from beam to beam. Currently, only 16 receivers are available at MAARSY. These 16 receive signals can be selected from groups of seven antennas, each called “hexagons”, or from a group of seven hexagons called “anemones” (see, e.g.,  Latteck et al., 2012b, for further technical details). For this campaign, we conducted an imaging experiment using 15 hexagons on reception similar to Sommer and Chau (2016)'s experiment. One receiver is always connected to the full antenna array and it is used as in the standard multiple experiments. The radar parameters of this experiment are summarized in Table 1.

Table 1 Parameters of MAARSY MIMO experiment.

2.2 MAARSY MIMO configuration

In order to improve the performance of our imaging experiment, we applied a coherent MIMO technique (Urco et al., 2018). The technique employs multiple independent transmitting antennas and multiple receiving antennas, both spatially separated, to take advantage of the transmit–receive geometry and to increase the angular resolution of the radar. If the antennas are closely separated or collocated, the signals from each transmitting–receiving path are coherent and can be combined to form a larger virtual receiving array. The resulting number of virtual receivers is equal to the number of transmitters times the number of receivers.

Depending on the transmitting and receiving antenna configuration, some virtual receivers can be redundant. In our experiment, we carefully selected the transmitting and receiving antenna configuration to get three special redundant virtual receivers. These three redundant virtual receivers were used for phase calibration of the transmitters as was done by Urco et al. (2018). Figure 1a shows the 15 hexagons used in reception and the three anemones used in transmission (B, D, F). Figure 1d shows the resulting virtual receiving antennas, whereby three of them are redundant and located at the origin.

Figure 1 MAARSY antenna configuration for SIMO (a, b, c) and MIMO (d, e, f). (a) 16 hexagons used in reception are shown in grey and three anemones used in transmission are colored. (b) Visibility samples for SIMO. (c) Point spread (or instrument) function for SIMO. (d) The virtual position of the resulting receiving antennas by using MIMO. (e) Visibility samples for MIMO. (f) Point spread (or instrument) function for MIMO. See text for further details.

In order to separate the contribution of each transmitter, a form of transmit diversity was needed. In Urco et al. (2018), three types of transmit diversity were proposed: code, time, and polarization. Code diversity is recommended for atmospheric observations given that this is not sensitive to the temporal correlation or polarization of the target of interest. Unfortunately, code diversity cannot be currently used in MAARSY. For targets for which the temporal correlation is less than the time separation between transmitters, time diversity can be applied. Given that PMSEs have a relative long correlation time (a few hundreds of milliseconds), we applied time diversity to enhance the spatiotemporal features of PMSE. The effective time separation between transmitters was 4, 4, and 8 ms between pairs BD, DF, and BF, respectively.

As explained by Urco et al. (2018), in a monostatic coherent MIMO radar, the relationship between the normalized spatial cross-correlation of signals from two different transmitting-receiving paths and the angular distribution of scattered power for a given range and frequency bin can be described by

$$\begin{array}{ll}\text{(1)}& {\displaystyle \frac{\langle v_{m,p}\cdot v_{n,q}^{*}\rangle }{\sqrt{\langle {\left|v_{m,p}\right|}^{\mathrm{2}}\rangle \langle {\left|v_{n,q}\right|}^{\mathrm{2}}\rangle }}}& {\displaystyle }=V\left(\mathit{k}\left(\mathrm{\Delta }{\mathit{r}}_{m,n}+\mathrm{\Delta }{\mathit{r}}_{p,q}\right)\right){\displaystyle }& {\displaystyle }=e^{j\left(\mathrm{2}\mathit{\pi }{f}_{\mathrm{d}}\mathit{\tau }+{\mathit{\varphi }}_{m,n}+{\mathit{\varphi }}_{p,q}\right)}\cdot \int B\left(\mathit{\theta }\right)\\ \text{(2)}& {\displaystyle }& {\displaystyle }\cdot e^{-j\mathit{k}\left(\mathrm{\Delta }{\mathit{r}}_{m,n}+\mathrm{\Delta }{\mathit{r}}_{p,q}\right)}\mathrm{d}\mathit{\theta },\end{array}$$

where v_m, p is the signal from the transmitting-receiving path (m,p), m being the receiver and p the transmitter; $v_{n,q}^{*}$ is the complex conjugate of the signal from the transmitting-receiving path (n,q), n being the receiver and q the transmitter; $\langle v_{m,p}*v_{n,q}^{*}\rangle$ is the cross-correlation of two signals from antennas spatially separated; $V\left(\mathit{k}\left(\mathrm{\Delta }{\mathit{r}}_{m,n}+\mathrm{\Delta }{\mathit{r}}_{p,q}\right)\right)$ is the visibility sample at $\left(\mathrm{\Delta }{\mathit{r}}_{m,n}+\mathrm{\Delta }{\mathit{r}}_{p,q}\right)$; k is the wave number vector equal to (2π∕λ)θ. And λ is the radar wavelength; θ is the angle of arrival equal to $({\mathit{\theta }}_x,{\mathit{\theta }}_y,{\mathit{\theta }}_z)$ which are the direction cosines in the (x, y, z) direction; B(θ) is the angular scattered power distribution, also known as brightness; Δr_m, n is the spatial separation between receivers m and n; Δr_p, q is the spatial separation between transmitters p and q; 2πf_dτ is the phase difference due to the Doppler shift of the target, f_d, where τ is the time separation between transmitters; ϕ_p, q is the phase difference between transmitters; ϕ_m, n is the phase difference between receivers.

A quick comparison between the visibility (sampling domain) for SIMO and MIMO, shown in Fig. 1b and e, indicates that the antenna aperture for MIMO is larger than the SIMO by ∼50 %. The difference lies in that the MIMO antenna aperture is defined as the maximum separation between two virtual receiving antennas, i.e., $\text{Max}\left(\mathrm{\Delta }{\mathit{r}}_{m,n}+\mathrm{\Delta }{\mathit{r}}_{p,q}\right)$, whereas for SIMO, $\mathrm{\Delta }{\mathit{r}}_{p,q}=\mathrm{0}$, and the antenna aperture is only defined by the maximum spatial separation between two receiving antennas. Figures 1c and f show the resulting instrument function or point spread function for SIMO and MIMO, respectively. As expected the half-power beamwidth (HPBW) for MIMO is ∼50 % smaller than for SIMO, resulting in an angular resolution of 2.4^∘ for MIMO compared to 3.6^∘ for SIMO. Furthermore, the sidelobes in the MIMO configuration are strongly reduced, given that the visibility is larger and contains no gaps.

3 Radar imaging implementation

Before inverting Eq. (1), the three phase differences due to time diversity (2πf_dτ), to receivers (ϕ_m, n), and to transmitters (ϕ_p, q) need to be corrected. When the analysis is done in the frequency domain we can easily correct the value 2πf_dτ, given that we know the frequency and the time separation between transmitters. On the other hand, the phase offsets between receivers have been calibrated using Cassiopeia A as a radio source (e.g., Chau et al., 2014). Additionally, we have calibrated the phase offset between transmitters using the three redundant virtual receivers described above. Each of the redundant virtual receivers comes from one transmitter. They were compared to have zero phase difference between each other, given that the three of them must be located in the same virtual position (see, e.g., Urco et al., 2018, for more details).

Once the imaging system is calibrated we can invert Eq. (1) to obtain the estimated brightness $\widehat{B\left(\mathit{\theta }\right)}$. Given that the number of unique visibility samples is still less than the number of unknowns (brightness points), some kind of regularization is needed to solve Eq. (1). Two of the most well-known radar imaging techniques applied to atmospheric and ionospheric targets are Capon (Palmer et al., 1998) and MaxEnt (Hysell, 1996).

3.1 Capon technique

As described by Kudeki and Sürücü (1991), the angular resolution obtained from a direct inversion of Eq. (1) using the inverse Fourier transform is limited by the longest baseline and the unmeasured antenna separations (visibility gaps). Palmer et al. (1998) proposed a new technique to improve the angular resolution based on the work of Capon (1969). Capon can be seen as an extension of the inverse Fourier transform. The difference lies in the fact that Capon chooses the antenna weights adaptively in order to minimize the sidelobe interference from signals outside of the direction of interest according to the data. Capon's technique provides an estimate of the brightness function given by

$$\begin{array}{}\text{(3)}& {\displaystyle }\widehat{B\left(\mathit{\theta }\right)}={\displaystyle \frac{\mathrm{1}}{M^{\mathrm{H}}\cdot V^{-\mathrm{1}}\cdot M}},\end{array}$$

where

$${\displaystyle }M={\left[e^{-jk\left(r_{m_{\mathrm{0}}}+r_{p_{\mathrm{0}}}\right)},e^{-jk\left(r_{m_{\mathrm{0}}}+r_{p_{\mathrm{1}}}\right)},\mathrm{\dots },e^{-jk\left(r_{m_{\mathrm{1}}}+r_{p_{\mathrm{0}}}\right)},\mathrm{\dots }\right]}^T$$

is the Fourier kernel and $V=V\left\{k\left(\mathrm{\Delta }{r}_{m_i,m_j}+\mathrm{\Delta }{r}_{p_k,p_l}\right)\right\}$ is the visibility due to the virtual receivers $v_{m_i,p_k}$ and $v_{m_j,p_l}$, with i and j being the receiver indices and l and k being the transmitter indices.

3.2 Maximum entropy technique

Even when MIMO is used, the problem is still underdetermined. Thus, there are infinite possible image solutions, B, which agree with the data, V. Of all possibilities, MaxEnt chooses the solution with the maximum entropy or minimal information content (e.g., Hysell, 1996) as the one to be the most likely brightness distribution and the most consistent with the available visibility data and their statistical uncertainties. The entropy for a given frequency bin and range can be defined as

$$\begin{array}{}\text{(4)}& {\displaystyle }& {\displaystyle }S=\sum _{\mathit{\theta }}\widehat{B\left(\mathit{\theta }\right)}\ln \widehat{\left\{B\left(\mathit{\theta }\right)/F\right\}},\text{(5)}& {\displaystyle }& {\displaystyle }F=\sum _{\mathit{\theta }}\widehat{B\left(\mathit{\theta }\right)},\end{array}$$

where F is the summation of the brightness distribution over the region of interest. The solution of Eq. (1) is defined by

$$\begin{array}{}\text{(6)}& {\displaystyle }\underset{\mathit{\theta }}{max}\left\{S\right\}\text{subject to}\left|V-M\cdot \widehat{B\left(\mathit{\theta }\right)}\right|<\mathit{ϵ},\end{array}$$

where ϵ is the noise amplitude associated with the visibility measurements. In this work, we have also considered the improvements of Hysell and Chau (2006). Specifically, we have taken into account the transmitting beam pattern and the statistical uncertainties of all the visibility pairs.

4 Results

Figure 2 shows the resulting 24-bit range–time Doppler intensity (RTDI) image of the vertical beam for 32 h of continuous operation on 16 and 17 July 2017. This plot was obtained after applying MaxEnt to the data and selecting the values that belong to the zenith angle. The signal intensity is represented as lightness, Doppler information as hue, and spectral width as saturation. As shown later, the resulting HPBW for this experiment is <1^∘, indicating that the Doppler information must be mainly due to the vertical motion. The RTDI plot indicates that the vertical motion is slow (green color) as expected. Nevertheless, there are two regions at 23:30 and 06:30 LT around 89 km in which the Doppler velocity presents unrealistic values. Indeed, PMSE were too strong at that time, so even the antenna sidelobes can be seen. Unfortunately, the imaging algorithm cannot assign the correct angle of arrival to these unusually strong echoes due to the angular ambiguity associated with our antenna array. The angular ambiguity is defined by the minimum separation between two antennas. The smaller the separation, the larger the angle without ambiguity (e.g., Woodman, 1997). A manual angular correction can be applied knowing the Doppler but it is a hard task in the presence of many targets. A smaller baseline is recommended in future experiments for these special cases.

Figure 2 A 24-bit image of a range–time Doppler intensity (RTDI) plot of PMSEs using MIMO, with time diversity conducted on 16 and 17 July 2017. The signal intensity is represented as lightness, Doppler information as hue, and spectral width as saturation. The legend on the left represents the SNR vs. Doppler color map for a saturation of 90 %. The legend on the right represents the spectral width vs. Doppler for a lightness of 50 %. Note that only the signal corresponding to the narrow region in the illuminated area is shown.

4.1 SIMO vs. MIMO results

Since the estimated brightness is expressed in polar coordinates $({\mathit{\theta }}_x,{\mathit{\theta }}_y,R)$, a cubic spline interpolation was applied to convert them to Cartesian coordinates, $B({\mathit{\theta }}_x,{\mathit{\theta }}_y,R)$ to $B(x,y,z)$, with the radar being located at the center (x=0, y=0, z=0). Below we show the results of two selected events (Events 1 and 2) after performing such interpolation. For both events, we show x vs. y cuts for a given z, as well as x vs. z cuts for a given y, where x, y, and z represent the east–west (EW) direction, north–south direction (NS), and altitude, respectively.

Examples of EW–NS and EW–altitude 2-D images for Event 1 obtained by applying Capon and MaxEnt to two different antenna configurations, SIMO and MIMO, are shown in Figs. 3 and 4, respectively. Four different results are shown: (a) SIMO-Capon, (b) SIMO-MaxEnt, (c) MIMO-Capon, and (d) MIMO-MaxEnt. Having a quick look at the results, it is clearly observable that (1) MaxEnt outperforms Capon when the same antenna configuration is used, either SIMO or MIMO. This was already pointed out by previous works (e.g., Yu et al., 2000; Harding and Milla, 2013). (2) As expected, MIMO shows a cleaner and more defined image compared to SIMO, when either Capon or MaxEnt is employed. (3) The improvement of using MIMO instead of SIMO is much better in MaxEnt than Capon. The improvement of MIMO-Capon with respect to SIMO-Capon is about 50 % due to the larger virtual antenna array, whereas the improvement of MIMO-MaxEnt with respect to SIMO-MaxEnt is much better than 50 %. This difference lies in the fact that Capon tries to reduce the sidelobes, adaptively steering them to echo-free zones. Unfortunately, for the two events shown, most of the illuminated area is filled with PMSE scattering, and thus the performance of Capon is expected to be comparable to the conventional beamforming (inverse Fourier transform). In the case of MaxEnt, the improvement is mainly due to the larger virtual antenna array and the use of statistical uncertainties as described by Hysell and Chau (2006). Unlike Capon, our MaxEnt implementation takes advantage of the redundant visibility pairs, giving more weight to pairs with less uncertainty, i.e. more redundancy. Figure 1b and e show the visibility pairs and their redundancy for SIMO and MIMO, respectively.

Figure 3 EW–NS images for z=85.8 km obtained from applying four different implementations, (a) SIMO-Capon, (b) SIMO-MaxEnt, (c) MIMO-Capon, and (d) MIMO-MaxEnt, at 00:56:55 UT on 17 July 2017, i.e., Event 1. Images are color coded the same as in Fig. 2. The yellow dashed horizontal and vertical lines represent the location of the NS–EW cuts shown in later figures for Event 1.

Figure 4 Similar to Fig. 3, but for an EW–altitude cut at y=0 km. The yellow dashed horizontal lines represent the location of the altitude cuts shown in previous and later figures for Event 1.

Coming back to our comparison of SIMO vs. MIMO, with MIMO-MaxEnt, small wave-like structures of 2 km wavelength can be clearly observed, which are invisible in SIMO implementations or MIMO-Capon. For example, observe the two wavefronts at x=10 km in Fig. 3d, right beside the larger meridionally oriented wavefronts of 7 km wavelength. This indicates that wave-like structures of different wavelengths coexist within PMSEs as previously seen in NLCs (e.g., Baumgarten and Fritts, 2014). In addition, Fig. 4d shows that the ascending structures (red color) have a higher signal-to-noise ratio (SNR) than the descending structures (blue color).

We show similar 2-D cuts for Event 2 in Figs. 5 and 6 for z=82.7 km and $y=-\mathrm{6}$ km, respectively. In this case, the observed wavelength is 12 km. Unlike the first event, the SNR is similar for targets with negative and positive Doppler. Figure 6d shows two very interesting points: (a) a very well-defined wave-like structure between 82 and 84 km and (b) a quasi-uniform structure between 84 and 86 km, which apparently has been modulated by the first wave. In this case, the wave-like structure is easily discernible even with SIMO-Capon, given that the wavelengths are larger than in Event 1 (see Fig. 5a).

Figure 5 Same as Fig. 3 but at 05:56:13 UT on 17 July 2017, i.e., Event 2. The yellow dashed horizontal and vertical lines represent the location of the NS/EW cuts shown in later figures for Event 2.

Figure 6 Same as Fig. 5 but for an EW–altitude cut at $y=-\mathrm{6}$ km. The yellow dashed horizontal lines represent the location of the altitude cuts shown in previous and later figures for Event 2.

4.2 MIMO results

Having shown the better qualitative performance of MIMO-MaxEnt with respect to the other three implementations for the two selected events above, next we present extended results using just MIMO-MaxEnt.

Figure 7 shows the evolution in time of two selected events, i.e., Event 1 (panels a, b, and c), and Event 2 (panels d, e, and f). Figure 7a and d show the time evolution vs. altitude for selected EW and NS coordinates. In these plots, we can appreciate how variable PMSE structures are, showing different altitudinal extensions. Note that the effective horizontal area is less than 1 km² in both cases.

Figure 7 24-bit time representation images of PMSE structures as a function of altitude (RTDI) (a, d), EW location (keogram) (b, e), and NS location (keogram) (c, f) for selected cuts, for both Event 1 (a, b, c) and Event 2 (d, e, f). In the keograms, the wind components obtained with specular meteor radars (SMRs) and MAARSY PMSEs are shown by pink and yellow arrows, respectively. The white dashed horizontal lines represent the location of the altitude, EW and NS cuts shown in previous figures, and current keograms for Events 1 and 2. The white dashed vertical lines represent the time of the cuts shown in previous figures.

The second and third row of Fig. 7 show the time evolution vs. EW direction and the time evolution vs. NS direction, for EW and NS keograms, respectively. We have included the zonal (u₀) and meridional (v₀) wind velocity estimated from combining a couple of specular meteor radars (SMRs) (pink arrow) and from MAARSY based on PMSE Doppler velocities (yellow arrow). The wind values are shown in Table 3. Since this is a “time” vs. “distance” plot, the zonal and meridional wind are represented by arrows for which their slopes indicate the wind magnitude, i.e, how long a target is displaced in the y axis for a certain time in the x axis. The SMR winds were obtained from combining SMR detections from Andenes and Tromsø in northern Norway (see, e.g., Chau et al., 2017, for details). In order to estimate the winds from PMSEs we used the following formula:

$$\begin{array}{}\text{(7)}& {\displaystyle }{v}_{\mathrm{rad}}\left({\mathit{\theta }}_x,{\mathit{\theta }}_y,{\mathit{\theta }}_z\right)=u\cdot {\mathit{\theta }}_x+v\cdot {\mathit{\theta }}_y+w\cdot {\mathit{\theta }}_z,\end{array}$$

where v_rad is the radial wind, $({\mathit{\theta }}_x,{\mathit{\theta }}_y,{\mathit{\theta }}_z)$ are the direction cosines, and (u, v, w) are the zonal, meridional, and vertical wind direction, respectively. Assuming a constant u, v, and w for a given altitude bin and time bin, and taking all the measurements with a SNR higher than −5 dB, we invert Eq. (7) and get u₀, v₀, and w₀, mean values of u, v, and w respectively.

The keograms for Event 1, i.e., from 00:50 to 01:05 UTC, show that the meridionally oriented wavefronts have a limited vertical extent centered at 85 km (Fig. 7a), and since this wave has a finite wavelength in the EW direction, the zonal wave propagation can clearly be observed in Fig. 7b, which shows that the elongated meridionally oriented wavefronts are zonally drifting with the same direction and with the same speed as the wind. In the NS direction, the meridional drifting of the wave is not clearly observed due to the elongated structure. Mesospheric wave-like features observed with airglow imagers (ripples) have also been noticed to drift with the background wind (e.g., Hecht, 2003). These ripples have been associated with gravity wave breaking and are a clear signature of atmospheric instability.

Figure 7d shows another interesting wave-like example (Event 2). Unlike the first case, this wave does not keep its amplitude in the vertical direction; see Fig. 7d. It grows and then disappears. Its direction of propagation in the zonal and meridional direction is also interesting. As shown in Fig. 7e, the direction of propagation in the zonal direction is completely opposite to the background wind. Whereas the wind is going from east to west, the wave propagates from west to east. In the NS direction (Fig. 7f), the wind is close to zero and we do not expect changes in this direction. Since its wavelength is relatively small, this structure might be classified as an instability; however, the opposite direction of propagation suggests that it could be a propagating gravity wave. Further investigation of these events is needed to understand the physical mechanisms behind them, including lidar and airglow imager observations.

PMSEs have been used as a neutral wind tracer, assuming that u, v, and w are constant and homogeneous during the analyzed time (e.g., Balsley and Riddle, 1984; Fritts et al., 1990; Hoppe and Fritts, 1995; Stober et al., 2013). Therefore, those works assumed that scatters from PMSEs are moving with the neutral wind at the same velocity and in the same direction. Unlike winds obtained from SMR, winds from PMSEs are affected by local disturbances, as shown in Fig. 7e and f. When the dynamics of local structures are not in agreement with the wind dynamics, a bias could be introduced in the wind estimation (as shown in the Event 2). However, when these local disturbances are moving with the wind, the estimated wind is not affected (Event 1). Note that the PMSE winds are in good agreement with the SMR winds in Event 1, but they are not for Event 2, particularly for the meridional component.

An animated sequence of the two events has been included in the Video supplement, i.e., Movies S1 and S2. For both events, the sequence includes selected cuts of EW–NS, EW–altitude, and NS–altitude. In Movie S1, we identify at least four examples of monochromatic waves with different wavelengths drifting with the wind in the northwest direction (at 23:57:37, 00:02:24, 00:10:57, and 00:55:33 UTC). Interestingly, in this case, longitudinal and transverse waves both drift with the background wind. In Movie S2, we show the complete evolution in time of Event 2. In the EW–altitude cut, the wave structure between 82 and 85 km drifts against the wind, whereas a layer at 87 km between 05:20 and 05:30 UTC follows the background wind. Note the projected radial wind (from red to blue) indicates a westward wind. These events are good examples of the complicated dynamics within PMSEs. Further analysis and interpretation of these high-resolution spatiotemporal structures will be done in a future work.

Figures 8 and 9 show 3-D maps of (a) the signal-to-noise ratio (SNR), (b) radial velocity, (c) locally enhanced SNR, and (d) residual radial velocity (i.e., v_res), for Events 1 and 2, respectively. In addition contours of locally enhanced SNR are overplotted on both the radial velocities. The SNR and radial velocity were obtained from the first and second spectral moments (e.g., Doviak and Zrnić, 1993). The locally enhanced SNR has been obtained using a 2-D Gaussian function kernel with a width of 6 pixels. The local enhancements allow us to observe weak structures within the strong one. For example, wave fronts are distinguishable in Fig. 8c which were not visible in Fig. 8a. On the other hand, the residual radial velocity was estimated by removing the contributions of the estimated mean horizontal velocities in the measured radial velocities; i.e.,

$$\begin{array}{}\text{(8)}& {\displaystyle }{v}_{\mathrm{res}}\left({\mathit{\theta }}_x,{\mathit{\theta }}_y,{\mathit{\theta }}_z\right)=v_{\mathrm{rad}}\left({\mathit{\theta }}_x,{\mathit{\theta }}_y,{\mathit{\theta }}_z\right)-\left(u_{\mathrm{0}}\cdot {\mathit{\theta }}_x+v_{\mathrm{0}}\cdot {\mathit{\theta }}_y\right).\end{array}$$

Figure 8 The 3-D contour plots at 00:55:33 UT on 17 July 2017, i.e., Event 1, for four selected altitudes: 84, 84.6, 85.2, and 85.8 km. The following are shown for each altitude: (a) SNR, (b) radial velocity, (c) locally enhanced SNR, and (d) residual radial velocity. Contours on locally enhanced SNR are overplotted in both velocity plots.

Figure 9 Same as Fig 8, but at 05:54:01 UT on 17 July 2017 for altitudes 82, 82.7, 83.4, and 84 km., i.e., Event 2.

Assuming that the v_res is mainly due to the vertical motion, we can clearly see in Fig. 8d how up (red) and down (blue) structures drift across the illuminated area, maybe due to KHI. Similarly, Fig. 9 shows animated images of Event 2. In this case, the horizontal wind was small and most of the radial velocity was due to the vertical motion; i.e., radial velocity and residual velocities are almost the same. As mentioned above, in this event, the waves propagate horizontally against the weak horizontal wind.

The animated versions of Figs. 8 and 9 are shown in Movies S3 and S4, respectively. Although the information might be redundant when compared to Movies S1 and S2, we have decided to include them to provide a more standard view of typical spectral parameters of a multibeam radar.

Making a quantitative comparison between SIMO and MIMO for real targets is not an easy task. We need a prior knowledge of the brightness to make a good analysis. This is not the case for PMSE. Fortunately, our observations include echoes from specular meteors; see the bright echoes located at (−10.5, −12.5) in Fig. 3d. Indeed, meteor echoes can be observed in the PMSE region, but the great majority of them occur outside this window. When a meteor echo occurs in the PMSE altitude it will be short-lived (less than a few hundred milliseconds). In previous studies meteor echoes were treated as outliers and were removed from the measurements (e.g., Hashimoto et al., 2014). For our benefit they can also be used to evaluate the angular resolution that can be achieved with our implementations quantitatively. A specular meteor echo could be considered to be a point target. Along with its trajectory, the trail is long (hundreds of meters to a few kilometers), but its angular response is narrow. In the transverse direction to the trail, it is very narrow and its angular response is also narrow.

In Fig. 10 we show the normalized angular scattered power distribution for a specular meteor using SIMO and MIMO in combination with Capon and MaxEnt. As expected, the range resolution does not change for SIMO or MIMO (see Fig. 10a). We see a peak at 89.1 km and low power at other ranges. However, when comparing Capon and MaxEnt, MaxEnt shows us a clean power distribution along all the ranges, while Capon shows us a remaining sidelobe contamination in other ranges, coming from other angles. This indicates that, even with MIMO, Capon does not suppress the sidelobes as well as MaxEnt. Figure 10b and c show us the angular power distribution for θ_x and θ_y, respectively, for which the points are the samples for a given angle and the continuous line is a fitted Gaussian function. Using the fitted function, we estimated the half-power beamwidth (HPBW) for each implementation. Table 2 summarizes the angular resolution and the improvement factor for each method compared to the theoretical angular resolution of the full array MAARSY radar. As we expected the improvement between SIMO and MIMO is about 1.5, given that we increased the antenna aperture for MIMO by ∼50 %. When combining MIMO and MaxEnt, surprisingly, we got an angular resolution of ∼0.6^∘, i.e., more than 5 times better than MAARSY's HPBW.

Figure 10 Normalized angular power distribution of a specular meteor echo as a function of (a) range, (b) EW angle (θ_x), and NS angle (θ_y). The results are shown for all four implementations, i.e., SIMO-Capon (blue), SIMO-MaxEnt (orange), MIMO-Capon (green), and MIMO-MaxEnt (red).

Table 2 Performance of imaging techniques.

Table 3 Mean wind values for the two events presented.

5 Discussion

We have shown qualitatively and quantitatively that radar imaging of PMSEs is significantly improved when using MIMO instead of SIMO configurations, by at least 50 %. Two different imaging methods have been applied, Capon and MaxEnt. As expected from previous works, MaxEnt images are better than Capon images; however, MaxEnt is computationally more demanding. Similarly, we found that the quality of MIMO-Capon is comparable to SIMO-MaxEnt.

Even though MIMO allows us to improve the point spread function, it is not perfect. We expect some artifacts due to the sidelobes which are −15 dB weaker than the main lobe; see Fig. 1f. When strong and weak echoes coexist in the same region, some artifacts might be confused as weak echoes. Although Capon and MaxEnt help to minimize the sidelobe contribution, we are being conservative by employing a relatively large SNR threshold (>5 dB), i.e., discarding weak echoes which might be contaminated by strong sidelobe echoes. By doing this we are increasing the statistical significance of our results which are persistent in time and space, as shown in the animations and the keograms.

The preliminary results using MIMO-MaxEnt are allowing us to observe PMSEs with unprecedented horizontal resolution (<1 km) compared to multibeam scanning experiments (Stober et al., 2013), and therefore the identification of structures with horizontal wavelengths less than 10 km (e.g., Event 1 above). For structures with wavelengths of the order of 15–20 km or so, the other imaging implementations, i.e., SIMO-Capon, SIMO-MaxEnt, and MIMO-Capon, are sufficiently good to characterize them. These new capabilities will allow KHIs and general GWs (not only monochromatic) to be better identified and characterized than previously done at polar mesospheric heights during the summer. Our proposed technique complements previous observations that have been performed at nighttime when the sky is clear using airglows and lidars (e.g., Smith, 2013; Hecht et al., 2000, 2007; Taylor et al., 2007).

We will leave the detailed analysis and interpretation of these events and other events observed with this new capability for a future effort. In the following paragraphs, we discuss the technical results and propose future improvements.

The improved resolution using MIMO results from the larger effective visibility aperture and the larger number of independent samples, as compared to a SIMO configuration, i.e., 125 m instead of 76 m and 475 samples instead of 163 samples, respectively. In addition, the MaxEnt approach allows an improvement at least a factor of 2 in angular resolution compared to Capon. The maximum number of horizontal blobs that could theoretically be estimated for each range, time, and “color” (i.e., frequency bin) would be 79 ($=\mathrm{475}/\mathrm{6}$), whereby each blob is characterized by a 2-D Gaussian function with six parameters (e.g., Chau and Woodman, 2001). Another reason for the better results using MIMO-MaxEnt is the number of redundant visibility measurements. Although they do not provide additional information in terms of degrees of freedom, the redundancy helps to reduce the statistical uncertainties of such visibility samples. Recall in our MIMO implementation that there are 1980 visibility samples (45×44), and only 475 are independent.

Despite the significant improvement, not everything is positive about applying MIMO. In the following paragraphs, we discussed the critical points of applying MIMO in terms of (a) power-aperture reduction and (b) computational demands and real-time applicability.

As indicated by Urco et al. (2018), in atmospheric radars MIMO is applicable to targets with a large RCS, since a reduction of power aperture is inherent to MIMO. In our particular application to PMSEs, the transmitter sections were 1∕7 of the total area, and therefore also 1∕7 of the total transmitter power, i.e., −17 dB transmitting signal, of usual experiments. In reception, 15 groups of seven antennas (hexagons) were used instead of the 433 available antennas. Moreover, given the time multiplexing, the number of coherent integrations was reduced and therefore the noise was increased, when compared to standard operations. In total, the sensitivity of our MIMO experiment is 27 dB lower. Looking at the PMSE RCS in Fig. 2 of Latteck and Strelnikova (2015), our MIMO observations are limited to PMSEs with RCS larger than 10⁻¹⁴ m⁻¹, i.e., approximately 40 % of the usual seasonal MAARSY PMSE observations.

MaxEnt is known to be computationally more demanding than Capon in SIMO applications (e.g., Yu et al., 2000). In the case of MIMO, the computational demands are significantly increased, given the larger number of effective receivers, i.e., 45 instead of 15. In terms of visibility pairs, the increase is from 210 to 1980. In the case of Capon, real-time processing is still possible with these increased numbers of samples; however, MaxEnt for both SIMO and MIMO is not applicable in a real-time application. For example, for 80 s of data using an i5 PC with 15 cores, the processing times are 20 min and 3 h for SIMO-MaxEnt and MIMO-MaxEnt, respectively. A future improvement to make MIMO-MaxEnt faster would be to use only one value of each redundant visibility sample, i.e., to only work with 475 independent samples instead of all 1980 measured visibility samples. Such a value could be obtained either from the average of all the values sampling the same visibility or by preselecting only one of them. After all, many of the independent samples are obtained with only one sample (green dots in Fig. 1e).

In general, a critical point for PMSE imaging is the drifting nature of the echoes. PMSE correlation times are relatively short, and under stationary conditions, one would require a few minutes of incoherent integration to reduce the statistical uncertainties of the visibility estimates. However, the structures to be imaged might move between 2 and 5 km in 60 s for typical mesospheric motions (40–80 m s⁻¹), either from drifting with the background wind (Event 1) or from wave propagation (Event 2). These drifting structures limit further the angular resolution that can be accomplished by any method since the resulting image will be significantly blurred for integration times of a few minutes.

To deal with the drifting nature of PMSEs, in future studies we will explore tracking techniques, i.e., make use of this information to improve the angular resolution (e.g., Vaswani and Zhan, 2016). Given the computational demands of MaxEnt in particular when combined with MIMO, we will also explore radar imaging with compressed sensing (CS) techniques (e.g., Donoho, 2006; Candes and Wakin, 2008). Harding and Milla (2013) applied CS to Jicamarca F-region irregularities, and show that CS produces results similar to MaxEnt. Our plan is to use MIMO-MaxEnt as a reference for other radar imaging techniques using SIMO, for example, CS in combination with tracking. Besides the computational demands, MIMO might not be applicable at other atmospheric radar sites, and therefore the exploration of other techniques using SIMO is required.

An additional improvement to the current observations would be the use of shorter pulses and therefore better range resolution, for example, 150 m. Further improvement in range could also be accomplished by applying range imaging (e.g., Palmer et al., 1998; Yu and Palmer, 2001), particularly in combination with the radar imaging implementations of this work, allowing angular resolutions less than 1^∘.

6 Conclusions

In this work, we have successfully implemented coherent MIMO with radar imaging at MAARSY to observe PMSEs with unprecedented angular resolution. The obtained resolution results from the combination of a larger effective aperture, a higher number of independent visibility samples resulting from MIMO, and improved angular resolution resulting from MaxEnt. Quantitatively, the maximum angular resolution accomplished is ∼0.6^∘, which is equivalent to having a 450 m diameter visibility aperture at 53.5 MHz and is a significant improvement to the MAARSY standard angular resolution of 3.6^∘.

The preliminary results with MIMO-MaxEnt allowed us to clearly identify structures slightly less than 1 km in diameter and wave-like structures with horizontal wavelengths less than 10 km, with a time resolution around 60 s. The identification of such structures, with varying degrees of intensity, suggests that one has to be careful about using PMSEs for estimating the background wind assuming horizontal homogeneity. Not only is the vertical wind not homogeneous, but also the brightness is not homogeneous horizontally.

Given the relatively long temporal correlation of PMSEs, i.e., a few minutes, larger integration of the noisy visibility in time would allow fewer statistical uncertainties in the resulting images of the two events presented. However, PMSE structures drift as they are imaged; therefore long integration times result in angular smearing. In the future, we plan to use the drifting information to improve the angular resolution by applying tracking techniques.

As mentioned above, the implementation of MIMO-MaxEnt is computationally intensive and is currently not applicable to real-time processing. On the other hand, MIMO-Capon can be implemented in real-time processing. Our strategy for near-future observations would be to use MIMO-Capon for real-time processing and use MIMO-MaxEnt for special events until more efficient implementations and/or faster computers are available.