A Hybrid Beamforming Design for Massive MIMO LEO Satellite Communications

1 Introduction

Integrated satellite-terrestrial cellular networks are currently being pursued for long awaited commercial applications. The main expected benefit is an uninterrupted coverage, even in unserved/underserved areas. LEO constellations seem to be the most promising platforms for satellite-based non-terrestrial networks (NTNs), due to their relatively shorter propagation delay (Kodheli et al., 2017). Although this propagation delay between ground terminals and LEO satellites (3 ms at 1,000 km above the Earth surface) is non-negligible, LEO constellations can support continuity and ubiquity of the radio services, and also latency critical applications with requirements within tens of ms.

Beam configuration and reconfiguration in a link between a LEO satellite and a ground mobile terminal is challenging due to the long round trip delay, fast movement of the satellite and limited on board processing capabilities. Current LEO constellations (del Portillo et al., 2019) make use of fixed analog beams that illuminate a given area of the Earth’s surface, without any capability to steer narrow beams in the user directions. Digital precoding stages are designed independently of the analog beam with the purpose of reducing inter-beam interference among the beams illuminating adjacent regions of interest (ROI) (Devillers and Pérez-Neira, 2011). Moreover, these solutions assume that a single satellite is illuminating a large specific ROI, and handover procedures are put in place to enable satellite switching at the user side before the satellite movement causes loss of coverage.

Designing the LEO satellite footprint and the specific beam codebook at the satellite side in an integrated satellite-terrestrial network is challenging. The potential beampatterns need to be designed to guarantee coverage of the Earth’s surface, limit the inter-beam interference, maximize system throughput, and, at the same time, exhibit some degree of compatibility with conventional codebooks defined in cellular standards. Commercial efforts on deploying LEO constellations, not necessarily for integration with the terrestrial network, are described in (del Portillo et al., 2019; Xia et al., 2019). Sizes of the beam footprints for the different deployments are provided, but the details of the beam codebook are not available. The massive MIMO LEO satellite communication system designed in (You et al., 2020) considers the set of all feasible beams generated by a uniform planar array at the satellite rather than a codebook, without trying to reduce the inter-beam interference or the amount of on-board processing.

In this paper, we propose and evaluate a massive MIMO LEO satellite communication system operating in the Ku band, based on a hybrid beamforming architecture. In particular, we design first the footprint associated to a given LEO satellite. Then, we design the beam codebook for the hybrid beamforming stage such that the derived footprint is covered by the beams in the codebook. We build our design using a 2-D DFT-based grid of beams as in Type I and Type II 5G New Radio CSI codebooks (Miao et al., 2018; Ericsson, 2017a; Ericsson, 2017b; Samsung, 2017; Technical Specification Group Radio Access Network, 2017a; b), considering an oversampling factor to shape the available beampatterns to the size of the satellite footprint. Assuming this codebook at the satellite and commercially available SatCom antennas for the mobile terminals, we evaluate the coverage and throughput of this LEO massive MIMO system operating in the Ku band. Additional performance metrics such as SNR degradation over time or inter-beam interference are also evaluated. Numerical results provide the baseline performance of a state-of-the art massive MIMO LEO satellite communication system and allow the identification of open research challenges.

2 System Model

We consider the downlink of a massive MIMO LEO satellite communication system operating in the Ku band, as in current commercial deployments (del Portillo et al., 2019), and also compliant with the 3GPP proposal for 5G broadband satellite services (Angeletti and De Gaudenzi, 2020). We denote as $B_{\text{DL}}$ the bandwidth allocated to the downlink. The system model we propose is general, but in the performance evaluation section we will focus on the Communications on the move (COOM) 5G use case. A given LEO satellite in the constellation can cover an elliptical region of interest (ROI), with semiradius $R_{\text{x}}$ and $R_{\text{y}}$, using a set of potential beams, as illustrated in Figure 1. Thus, the system supports the direct transmission to a maximum number of simultaneous mobile user terminals (UTs) on the ground, denoted as $N_{\text{u}}$, by using $N_{\text{b}}$ analog beams, with $N_{\text{u}}<N_{\text{b}}$. Time or frequency division multiplexing mechanisms will enable beam sharing when more than one user per beam is simultaneously served. The beam footprints are not necessarily identical, to account for the different elevation angles across beams. We assume that the beams generated at the LEO satellite belong to a fixed codebook, i.e. the beams move with the satellite, as considered in Scenario C1 and D2 in the 3GPP technical report that specifies solutions for 5G NR to support NTNs (Technical Specification Group Radio Access Network, 2019). We also assume that the UTs know the position and trajectory of the satellites in the constellation at all times, as proposed in (Kim et al., 2020). This information is periodically updated and broadcasted by the LEO satellite operations center to the cellular network. Knowledge of the planned trajectory allows the UTs to obtain an estimate of the satellite position even at locations where the terrestrial network is not available. It is out of the scope of this paper to consider the impact of the position or trajectory errors in the system model.

FIGURE 1

FIGURE 1. System model for a LEO satellite communication system covering an elliptical area using a predefined beam codebook.

2.1 MIMO Architecture and Antenna Models

We assume that the satellite is equipped with a large uniform planar array (UPA) with dimension $N_{\text{x}}^{\text{sat}}\times N_{\text{y}}^{\text{sat}}$ and $N_{\text{x}}^{\text{sat}}\times N_{\text{y}}^{\text{sat}}$ denoting the total number of antennas. We consider $\text{x}$ to be the axis in the direction of the satellite’s movement and $\text{y}$ to be the axis in the direction orthogonal to the movement. The possible beamforming architectures for the LEO satellite include: fully digital, reconfigurable analog, analog based on beam switching and hybrid (Angeletti and De Gaudenzi, 2020). Our choice is a partially connected hybrid MIMO (Multi-Input Multi-Output) architecture with a fixed codebook-based analog stage, as a practical trade-off between complexity and reconfigurability. In this regard, a fully analogue reconfigurable beamforming network entails a large mass and power consumption, whereas a fully digital solution is costly and power hungry due to all the required connections. The advantages of the hybrid approach become especially relevant for very large arrays, with major efforts currently underway to develop massive hybrid arrays based on multiple subarrays operating with analog beam-forming (De Gaudenzi et al., 2020). The number of beamforming ports in each dimension is defined as $N_{\text{x}}^{\text{RF}}$ and $N_{\text{y}}^{\text{RF}}$ making a total of $N^{\text{RF}}=N_{\text{x}}^{\text{RF}}{N}_{\text{y}}^{\text{RF}}$ RF-chains. The large array is organized into sub-arrays of size $N_{\text{x}}^{\text{sub}}\times N_{\text{y}}^{\text{sub}}$, with $N_{\text{x}}^{\text{sub}}=N_{\text{x}}^{\text{sat}}/N_{\text{x}}^{\text{RF}},N_{\text{y}}^{\text{sub}}=N_{\text{y}}^{\text{sat}}/N_{\text{y}}^{\text{RF}}$ as illustrated in Figure 2. Each radio frequency (RF) chain drives one of the subarrays, generating an independent spot beam so the whole coverage area is illuminated. This approach allows the reduction of the number of beamforming ports and on board processing requirements, since the analog beams do not have to be constantly reconfigured.

FIGURE 2

FIGURE 2. Partially connected sub-array structure at the satellite antenna. Each subarray is controlled by a single RF-chain.

Assuming a hybrid architecture, the precoding matrix consisting on $N_{\text{b}}$ beam-patterns $\mathbf{F}\in {\mathrm{ℂ}}^{N^{\text{sat}}\times N_{\text{b}}}$ can be written as $\mathbf{F}={\mathbf{F}}_{\text{RF}}{\mathbf{F}}_{\text{BB}}$, where ${\mathbf{F}}_{\text{BB}}\in {\mathrm{ℂ}}^{N^{\text{RF}}\times N_{\text{b}}}$ is the digital precoder and ${\mathbf{F}}_{\text{RF}}\in {\mathrm{ℂ}}^{N^{\text{sat}}\times N^{\text{RF}}}$ is the analog precoder. Due to power constraints in the mth RF-chain, we have that ${\left|\left|{\left[{\mathbf{F}}_{\text{BB}}\right]}_{m,:}\right|\right|}^2\le P_{\text{RF}}$. We define the index function $I_{\text{RF}}:\left\{1,\dots ,N_{\text{sat}}\right\}\to \left\{1,\dots ,N^{\text{RX}}\right\}$, where $I_{\text{RF}}\left(n\right)$ is the index of the RF-chain to which the nth antenna is connected. With this definition, the RF precoder ${\mathbf{F}}_{\text{RF}}$ can be written as ${\left[{\mathbf{F}}_{\text{RF}}\right]}_{n,m}=\{\begin{array}{cc}{e}^{{\varphi }_{n}i}& \text{if }m=I_{\text{RF}}\left(n\right)\\ 0& \text{otherwise}\end{array}$. We envision a digital precoder at the satellite to manage inter-beam interference in the context of a full frequency reuse system. Finally, it is also assumed that all the power amplifiers at the satellite UPA are identical, and their impact can be modeled as a scaling factor $\sqrt{P}$ (Angeletti and De Gaudenzi, 2020).

For the antenna arrays at the mobile terminal, we consider three possible solutions that enable beamforming in the satellite direction. The first one is a uniform planar array (UPA) of moderate size, with dimensions $N_{\text{x}}^{\text{UT}}\times N_{\text{y}}^{\text{UT}}$, having a total of $N^{\text{UT}}=N_{\text{x}}^{\text{UT}}{N}_{\text{y}}^{\text{UT}}$ antenna elements. Recent advances on antenna technology enable the commercial applications of these planar arrays for SatCom terminals (Sikri and Jayasuriya, 2019). The second considered solution is a flat panel antenna based on leaky wave technology (Mehdipour and Eleftheriades, 2014). The final option is a commercial liquid crystal based metasurface antenna developed by Kymeta (Mehdipour et al., 2019).

We define $\mathbf{w}\in {\mathrm{ℂ}}^{N^{\text{UT}}}$ as the vector containing the UT antenna weights over the set of feasible configurations. We consider first the UPA with antenna elements relative positions defined through ${\mathbf{k}}_{\text{u}}\in {\mathrm{ℝ}}^{N^{\text{RF}}\times 3}$, being ${\left[{\mathbf{k}}_{\text{u}}\right]}_{n,:}$ the nth antenna element’s relative position. In this case, the UT steering vector ${\mathrm{a}}_{\text{u}}:{\mathcal{S}}_2\to {\mathrm{ℂ}}^{N^{\text{UT}}}$ is defined over the unitary sphere ${\mathcal{S}}_2$ as ${\left[{\mathbf{a}}_{\text{u}}\left(\mathbf{v}\right)\right]}_n=e^{-i\frac{2\pi }{\lambda }\langle {\left[{\mathbf{k}}_{\text{u}}\right]}_{n,:},\mathbf{v}\rangle }$. The maximum gain in a direction $\mathbf{v}\in {\mathcal{S}}_2$ can be computed as ${\text{max}}_{\mathbf{w}}\frac{\left|\langle {\mathbf{a}}_{\text{u}}\left(\mathbf{v}\right),\mathbf{w}\rangle \right|}{\left|\left|\mathbf{w}\right|\right|}$. The set of feasible beam-patterns for a uniform planar array is ${\left\{z\in \mathrm{ℂ}:\left|z\right|=1\right\}}^{N^{\text{UT}}}$, and its maximum achievable gain in any direction is $10{\text{log}}_{10}\left(N^{\text{UT}}\right)$ when each one of the antenna elements is assumed to provide a unitary gain.

For a leaky wave antenna with leakage factor α, a feasible beam-patterns $\mathbf{w}\in {\mathrm{ℂ}}^{N^{\text{UT}}}$ follows the expression ${\left[\mathbf{w}\right]}_n=e^{-\frac{1}{\lambda }\left({\left[{\mathbf{k}}_{\text{u}}\right]}_{n,1}\alpha +{\left[{\mathbf{k}}_{\text{u}}\right]}_{n,2}\alpha \right)-\frac{i}{\lambda }\left({\left[{\mathbf{k}}_{\text{u}}\right]}_{n,1}{\theta }_1+{\left[{\mathbf{k}}_{\text{u}}\right]}_{n,2}{\theta }_2\right)}$ for ${\theta }_1,{\theta }_2\in \left[-\pi ,\pi \right]$ (Mehdipour and Eleftheriades, 2014). When $\alpha \ne 0$, the gain of the leaky wave antenna can be computed as $10{\text{log}}_{10}\left(\frac{\left(1-e^{-\alpha N_{\text{x}}^{\text{UT}}}\right)\left(1-e^{-\alpha N_{\text{y}}^{\text{UT}}}\right)\left(1+e^{-\alpha }\right)\left(1+e^{-\alpha }\right)}{\left(1+e^{-\alpha N_{\text{x}}^{\text{UT}}}\right)\left(1+e^{-\alpha N_{\text{y}}^{\text{UT}}}\right)\left(1-e^{-\alpha }\right)\left(1-e^{-\alpha }\right)}\right)$. Note that there is a gain threshold $10{\text{log}}_{10}\left(\frac{\left(1+e^{-\alpha }\right)\left(1+e^{-\alpha }\right)}{\left(1-e^{-\alpha }\right)\left(1-e^{-\alpha }\right)}\right)$ that cannot be exceeded by increasing the number of antenna elements in the leaky-wave antenna.

Finally, the metasurface antenna designed by Kymeta has a peak gain of 33 dB and a 1.2 cosine roll-off factor. This translates into a gain of $33-1.2\cdot 10{\text{log}}_{10}\left(\text{sin}\left({\theta }_{\text{el}}\right)\right)$ dB towards a satellite with an elevation angle ${\theta }_{\text{el}}$ relative to the user. The analytical expression of the array factor for this antenna is not available, so only system metrics that depend on the peak gain can be computed when the terminal is operating with this antenna.

2.2 Channel Model

We consider the channel between a single LEO satellite located at position ${\mathbf{x}}_{\text{sat}}\in {\mathrm{ℝ}}^3$ and a UT located at position ${\mathbf{x}}_{\text{u}}\in {\mathrm{ℝ}}^3$ on the Earth’s surface, both measured from Earth’s center. The satellite is orbiting Earth at a height $h_{\text{sat}}$, moving with a linear speed of

$v_{\text{sat}}=\sqrt{G\frac{m_{\text{earth}}}{R_{\text{earth}}+h_{\text{sat}}}},(1)$

where $R_{\text{earth}}$ is the Earth’s radius, $m_{\text{earth}}$ is the Earth’s mass, and G is the gravitational constant. This is equivalent to an angular speed of $w_{\text{sat}}=\frac{v_{\text{sat}}}{r_{\text{earth}}+h_{\text{sat}}}$. The closest point to the satellite on the Earth surface, from now on the satellite’s shadow, moves at a linear speed of $v_{\text{shadow}}=w_{\text{sat}}{r}_{\text{earth}}$ and circles the Earth with a period of $\frac{2\pi }{w_{\text{sat}}}$.

2.2.1 Link Budget Model

The link power loss includes the free space path loss, denoted as $L{P}_{\text{fs}}$, and the atmospheric absorption loss $L{P}_{\text{at}}$. We neglect the impact of the atmospheric fading on the link budget. This way, the signal-to-noise ratio (SNR) can be written in terms of the received signal strength (RSS) and the noise power ${\sigma }^2$ as (Technical Specification Group Radio Access Network, 2019)

$SNR\left[\text{dB}\right]=RSS\left[\text{dBW}\right]-{\sigma }^2\left[\text{dBW}\right](2)$

$RSS\left[\text{dBW}\right]=P_{\text{TX}}\left[\text{dBW}\right]-L{P}_{\text{cable}}\left[\text{dB}\right]+G_{\text{TX}}\left[\text{dB}\right]-L{P}_{\text{at}}\left[\text{dB}\right]-L{P}_{\text{fs}}\left[\text{dB}\right]+G_{\text{RX}}\left[\text{dB}\right](3)$

${\sigma }^2\left[\text{dBW}\right]=T\left[\text{dBK}\right]+k[\text{dBW}/\text{K}/\text{Hz}]+B\left[\text{dBHz}\right],(4)$

where $P_{\text{TX}}$ is the transmit power, $G_{\text{TX}}$ is the transmit antenna gain, $L{P}_{\text{cable}}$ is the cable loss between the antenna and the transmitter, $G_{\text{RX}}$ is the receiver antenna gain, T is the noise temperature, k is the Boltzmann constant equal to $-228.6\text{dBW}/\text{K}/\text{Hz}$, and B is the channel bandwidth.

The free space path loss in dB is given by

$L{P}_{fs}=20\left({\text{log}}_{10}\left(\left|\left|{\mathbf{x}}_{\text{sat}}-{\mathbf{x}}_{\text{u}}\right|\right|\right)+{\text{log}}_{10}\left(f\right)+{\text{log}}_{10}\left(\frac{4\pi }{c}\right)\right),(5)$

with $\left|\left|{\mathbf{x}}_{\text{sat}}-{\mathbf{x}}_{\text{u}}\right|\right|$ being the traveled distance, f the frequency and c the speed of light.

Defining $A\left(f,T\left(h\right),P\left(h\right),\rho \left(h\right)\right)$ as the loss per meter for temperature $T\left(h\right)$, pressure $P\left(h\right)$ and humidity $\rho \left(h\right)$ at height h, the atmospheric path loss due to absorption is

$L{P}_{at}={\displaystyle {\int }_{h_{\text{min}}}^{h_{\text{max}}}A\left(f,T\left(h\right),P\left(h\right),\rho \left(h\right)\right)dh.}(6)$

The integral is performed over the efficient limits of the atmosphere, that is between the user’s height $h_{\text{min}}$ ($h_{\text{min}}=0$ for sea level) and the atmosphere thickness $h_{\text{max}}=81km$. The atmospheric temperature and pressure as a function of the height can be extracted using the 1976 Standard Atmospheric model depicted in Figure 3, while the corresponding loss per meter can be computed using tabulated functions implemented in available software packages like gaspl in MATLAB. The humidity density vertical profile $\rho \left(h\right)$ can significantly vary depending on meteorological conditions, and its physical properties (condensation and freezing) can barely exist past its freezing point ($<-0$°C) without the formation of clouds. Therefore, we will assume that there are no clouds and that the water vapor is $7.5\text{g}/{\text{m}}^3$ at sea level $h=0$ and proportional to the air density in the first 2.3 km of the atmosphere, where the temperature is still over the $0\text{deg}$.

FIGURE 3

FIGURE 3. 1976 Standard Atmospheric model showing temperature and pressure vertical profiles.

2.2.2 Satellite Angular Speed and Doppler Effects

The satellite movement at a speed ${\mathbf{v}}_{\text{sat}}\in {\mathrm{ℝ}}^3$ introduces a relative angular speed $w_{\text{rel}}$ with respect to the user and, for a signal at frequency f, a Doppler effect $\text{Δ}{f}_{\text{sat}}$ with expressions

$\begin{array}{cc}\text{Δ}{f}_{\text{sat}}=& <{\mathbf{v}}_{\text{sat}},\frac{{\mathbf{x}}_{\text{u}}-{\mathbf{x}}_{\text{sat}}}{\left|\left|{\mathbf{x}}_{\text{u}}-{\mathbf{x}}_{\text{sat}}\right|\right|}>\frac{f}{c},\\ w_{\text{rel}}=& \sqrt{1-|<\frac{{\mathbf{v}}_{\text{sat}}}{{\mathit{v}}_{\text{sat}}},\frac{{\mathbf{x}}_{\text{u}}-{\mathbf{x}}_{\text{sat}}}{\Vert {\mathbf{x}}_{\text{u}}-{\mathbf{x}}_{\text{sat}}\Vert }>{|}^2}\frac{{\mathit{v}}_{\text{sat}}}{\left|\left|{\mathbf{x}}_{\text{u}}-{\mathbf{x}}_{\text{sat}}\right|\right|}.\end{array}(7)$

Note that the angular speed is computed in absolute terms and it is not directly related to a parametrization of choice, like for example polar coordinates from the UT point of view. Therefore, it does not translate to variational speed of azimuth and elevation angles. We discourage the use of satellite tracking strategies using a polar parametrization to avoid sharp changes in the azimuth direction. For example, considering a unitary sphere centered in the origin in an arbitrary Cartesian coordinate system and its corresponding polar parametrization, a small change of direction from $[\mathit{ϵ},0,\sqrt{1-{\mathit{ϵ}}^2}]$ to $[-\mathit{ϵ},0,\sqrt{1-{\mathit{ϵ}}^2}]$ for any $\mathit{ϵ}\in \left[\mathrm{0,1}\right]$ as small as desired, involves a ${180}^{\circ }$ change in the azimuth angle, while the effective change of direction is only of $\text{arccos}\left(1-2{\mathit{ϵ}}^2\right)\left[\text{rad}\right]$, that for small values of ϵ is close to $2\mathit{ϵ}\left[\text{rad}\right]$.

The Doppler effect caused by the satellite’s movement is increasing and anti-symmetric in the movement’s direction, thus its maximum is located at the elliptical footprint boundary. By creating orbits with many satellites, we can reduce the maximum absolute Doppler effect, since the distance between satellites decreases, thus reducing the size of the covered area in the direction of movement. As for the relative angular speed, the effect is symmetric concave with a maximum in the satellite’s shadow. Therefore, the maximum is not affected by the satellite’s density. Using basic geometry, these maximum values are given by the following formulas

$\begin{array}{cc}\text{max}\left|\text{Δ}{f}_{\text{sat}}\right|=& \frac{\frac{R_{\text{earth}}{R}_{\text{x}}}{\sqrt{R_{\text{earth}}^2+R_{\text{x}}^2}}\sqrt{\frac{G{m}_{\text{earth}}}{R_{\text{earth}}+h_{\text{sat}}}}}{\sqrt{R_{\text{earth}}^2+{\left(R_{\text{earth}}+h_{\text{sat}}\right)}^2-2\frac{R_{\text{earth}}^2\left(R_{\text{earth}}+h_{\text{sat}}\right)}{\sqrt{R_{\text{earth}}^2+R_{\text{x}}^2}}}}\frac{f}{c},\\ \text{max}{w}_{\text{rel}}=& \frac{1}{h_{\text{sat}}}\sqrt{G\frac{m_{\text{earth}}}{R_{\text{earth}}+h_{\text{sat}}}},\end{array}(8)$

The maximum Doppler caused by the satellite’s movement can be simplified when $R_{\text{x}}\ll R_{\text{earth}}$ to

$\text{max}\left|\text{Δ}{f}_{\text{sat}}\right|=\frac{R_{\text{x}}}{h_{\text{sat}}}\sqrt{\frac{G{m}_{\text{earth}}}{R_{\text{earth}}+h_{\text{sat}}}}\frac{f}{c}=\frac{R_{\text{x}}{v}_{\text{sat}}}{h_{\text{sat}}}\frac{f}{c}.(9)$

The Doppler effect caused by the UT moving with the speed vector ${\mathbf{v}}_{\text{u}}\in {\mathrm{ℝ}}^3$, for a random $\mathbf{v}\in {\mathcal{S}}_2$ which we are assuming to be uniformly distributed, is given by $\text{Δ}{f}_{\text{u}}=\langle {\mathbf{v}}_{\text{u}},\mathbf{v}\rangle \frac{f}{c}$. The maximum average value for this expression is given by $\text{max}\left|\text{Δ}{f}_{\text{u}}\right|=v_{\text{u}}\frac{f}{c}$, while the average absolute value can be computed as $\mathbb{E}\left(\left|\text{Δ}{f}_{\text{u}}\right|\right)=\frac{v_{\text{u}}}{6}\frac{f}{c}$, being $v_{\text{u}}=\Vert {\mathbf{v}}_{\text{u}}\Vert$ the absolute UT speed. Analagously, the Earth rotation induces a movement in the UT to the East at a speed of $v_{\text{rot}}=w_{\text{rot}}{r}_{\text{earth}}\text{cos}\left(\varphi \right)$, or equivalently ${\mathbf{v}}_{\text{rot}}=w_{\text{rot}}\left[{\left[{\mathbf{x}}_{\text{u}}\right]}_2,-{\left[{\mathbf{x}}_{\text{u}}\right]}_1,0\right]$, for $w_{\text{rot}}=\frac{2\pi }{86164.0905}[\text{rad}/\text{s}]$ Earth’s angular speed and ϕ UT’s latitude coordinate. Thus, the Doppler effect induced by the UT movement due Earth’s rotation is given by $\text{Δ}{f}_{\text{rot}}=<{\mathbf{v}}_{\text{rot}},\frac{{\mathrm{x}}_{\text{sat}}-{\mathrm{x}}_{\text{u}}}{\left|\left|{\mathrm{x}}_{\text{sat}}-{\mathrm{x}}_{\text{u}}\right|\right|}>\frac{f}{c}$, this expression simplifies to $\text{Δ}{f}_{\text{rot}}=\frac{w_{\text{rot}}}{\Vert {\mathbf{x}}_{\text{sat}}-{\mathbf{x}}_{\text{u}}\Vert }\left({\left[{\mathbf{x}}_{\text{u}}\right]}_2{\left[{\mathbf{x}}_{\text{sat}}\right]}_1-{\left[{\mathbf{x}}_{\text{u}}\right]}_1{\left[{\mathbf{x}}_{\text{sat}}\right]}_2\right)\frac{f}{c}$ and is bounded in absolute value by $\frac{\frac{R_{\text{earth}}\left(R_{\text{earth}}+h_{\text{sat}}\right)r}{\sqrt{R_{\text{earth}}^2+r^2}}{w}_{\text{sat}}}{\sqrt{R_{\text{earth}}^2+{\left(R_{\text{earth}}+h_{\text{sat}}\right)}^2-2\frac{R_{\text{earth}}^2\left(R_{\text{earth}}+h_{\text{sat}}\right)}{\sqrt{R_{\text{earth}}^2+r^2}}}}\frac{f}{c}$ for r being the distance between the UT and the satellite’s shadow, with equality only when both UT and satellite are located on the equator line.

2.2.3 Channel Matrix

Due to the large distance and lack of medium interaction, a Rician geometric channel model with power β and Rician factor $K_{\text{r}}$ can be assumed (You et al., 2020). We will also assume that the Doppler effect can be corrected and removed from the channel matrix. The UPA at the satellite can be described by steering vectors in the direction $\mathbf{v}\in {\mathcal{S}}_2$, denoted as ${\mathbf{a}}_{\text{sat}}:{\mathcal{S}}_2\to {\mathrm{ℂ}}^{N^{\text{sat}}}$. The antenna element positions are defined as ${\mathbf{k}}_{\text{sat}}\in {\mathrm{ℝ}}^{N^{\text{sat}}\times 3}$ with ${\left[{\mathbf{k}}_{\text{sat}}\right]}_{n,:}$ being the satellite’s nth antenna element relative position. With this definition, the steering vectors can be written as ${\left[{\mathbf{a}}_{\text{sat}}\left(\mathbf{v}\right)\right]}_n=e^{-i\frac{2\pi }{\lambda }\langle \left[{\mathbf{k}}_{\text{sat}}\right],n,:,\mathrm{v}\rangle }$. An alternative notation considering the spherical coordinates for the azimuth (ϕ) and elevation (θ) angles can be obtained by substituting $\mathbf{v}=\left[\text{cos}\left(\varphi \right)\text{cos}\left(\theta \right),\text{sin}\left(\varphi \right)\text{cos}\left(\theta \right),\text{sin}\left(\theta \right)\right]$. This way, the steering vectors are written as $\left[{\mathbf{a}}_{\text{sat}}{\left(\varphi ,\theta \right)}_n\right]=e^{-i\frac{2\pi }{\lambda }{\left[{\mathbf{k}}_{\text{sat}}\right]}_{n,1}\text{cos}\left(\varphi \right)\text{cos}\left(\pi \right)+{\left[{\mathbf{k}}_{\text{sat}}\right]}_{n,2}\text{sin}\left(\varphi \right)\text{cos}\left(\theta \right)}$. Given the steering vectors, the equivalent channel matrix after Doppler correction is given by

$\mathbf{H}=\left(\sqrt{\frac{K_{\text{r}}\beta }{K_{\text{r}}+1}}{\mathrm{a}}_{\text{u}}\left(\frac{{\mathbf{x}}_{\text{sat}}-{\mathbf{x}}_{\text{u}}}{\Vert {\mathbf{x}}_{\text{sat}}-{\mathbf{x}}_{\text{u}}\Vert }\right)+\sqrt{\frac{\beta }{K_{\text{r}}+1}}{\mathbf{a}}_{\text{R}}\right){\mathbf{a}}_{\text{sat}}^{\text{H}}\left(\frac{{\mathbf{x}}_{\text{u}}-{\mathbf{x}}_{\text{sat}}}{\Vert {\mathbf{x}}_{\text{u}}-{\mathbf{x}}_{\text{sat}}\Vert }\right),(10)$

where ${\mathbf{a}}_{\text{R}}\sim \mathcal{CN}\left(0,\text{Σ}\right)$ is the Rician component satisfying $\text{trace}\left(\text{Σ}\right)={\Vert {\mathbf{a}}_{\text{u}}\left(\frac{{\mathbf{x}}_{\text{u}}-{\mathbf{x}}_{\text{sat}}}{\left|\left|{\mathbf{x}}_{\text{u}}-{\mathbf{x}}_{\text{sat}}\right|\right|}\right)\Vert }^2$. Defining $\gamma =\sqrt{\frac{K_{\text{r}}\beta }{K_{\text{r}}+1}}$ as the complex gain comprising all path loss effects and phase shift corresponding to the LoS component, the channel matrix can be written in a more compact way as

$\mathbf{H}=\gamma \left({\mathbf{a}}_{\text{u}}\left(\frac{{\mathbf{x}}_{\text{sat}}-{\mathbf{x}}_{\text{u}}}{\Vert {\mathbf{x}}_{\text{sat}}-{\mathbf{x}}_{\text{u}}\Vert }\right)+\frac{1}{\sqrt{K_{\text{r}}}}{\mathbf{a}}_{\text{R}}\right){\mathbf{a}}_{\text{sat}}^{\text{H}}\left(\frac{{\mathbf{x}}_{\text{u}}-{\mathbf{x}}_{\text{sat}}}{\Vert {\mathbf{x}}_{\text{u}}-{\mathbf{x}}_{\text{sat}}\Vert }\right).(11)$

3 System Design

3.1 Design of the Coverage Area

We consider a LEO constellation where satellites are distributed around $N_{\text{p}}$ orbital planes with an inclination ${\theta }_{\text{op}}$, each one with $N_{\text{s}}$ satellites. Each satellite illuminates an elliptical area as shown in Figure 1. From a geometric perspective, to parametrize the area covered by the satellite, we are choosing a stereographic projection respect to the Earth’s center as illustrated in Figure 4A. This allows us to use Cartesian coordinates while still considering Earth’s intrinsic geometry. Mathematically, this parametrization of the covered area translates the Cartesian point $\left(x,y\right)$ into the sphere point $R_{\text{earth}}\frac{R_{\text{earth}}{\mathbf{u}}_{\text{sat}}+x{\mathbf{u}}_{\text{x}}+y{\mathbf{u}}_{\text{y}}}{\sqrt{R_{\text{earth}}^2+x^2+y^2}}$ for the unitary vector definition ${\mathbf{u}}_{\text{sat}}=\frac{{\mathbf{x}}_{\text{sat}}}{\left|\left|{\mathbf{x}}_{\text{sat}}\right|\right|}$, ${\mathbf{u}}_{\text{x}}=\frac{{\mathbf{v}}_{\text{sat}}}{\left|\left|{\mathbf{v}}_{\text{sat}}\right|\right|}$ and ${\mathbf{u}}_{\text{y}}={\mathbf{u}}_{\text{x}}\wedge {\mathbf{u}}_{\text{sat}}$ using $\wedge$ to denote the cross product.

FIGURE 4

FIGURE 4. Design of the satellite elliptical footprint. (A) Stereographic parametrization of the area covered by the satellite. (B) Illustration of the design process.

To design the sizes of the semiradiuses of the elliptical ROI, we will force the LEO constellation to cover the entire Earth. Thus, the equivalent flat problem we want to solve is to cover the entire Earth surface with ellipses centered in the vertices of a rectangular grid with spacing $\pi R_{\text{earth}}/N_{\text{s}}$ in the first dimension and $\pi R_{\text{earth}}\text{sin}\left({\theta }_{\text{op}}\right)/N_{\text{p}}$ in the second dimension. Due to ellipses being convex and symmetric in both axes, this is equivalent to design ellipses containing the center points of the rectangular grid, as illustrated in Figure 4B. This condition can be written as

${\left(\frac{\pi R_{\text{earth}}/N_{\text{s}}}{R_{\text{x}}}\right)}^2+{\left(\frac{\pi R_{\text{earth}}\text{sin}\left({\theta }_{\text{op}}\right)/N_{\text{p}}}{R_{\text{y}}}\right)}^2\le 1.(12)$

An example of a semiradius pair that fulfills this condition is

$R_{\text{x}}=\sqrt{2}\pi R_{\text{earth}}/N_{\text{s}},R_{\text{y}}=\sqrt{2}\pi R_{\text{earth}}\text{sin}\left({\theta }_{\text{op}}\right)/N_{\text{p}}.(13)$

3.2 Beam Codebook Design

Our goal is to design a set of analog beampatterns or beam codebook, such that the set of beam footprints covers an ellipsoidal ROI as that described in the previous section. Depending on the location of the active users, some active beams will be selected from the set available in the codebook. We will consider a DFT-based analog precoder codebook structure, a popular solution when operating with planar arrays. In this case, the array response vectors in each dimension have a linearly increasing phase, like the columns in the DFT matrix. Therefore, the MIMO channel can be written in terms of array response vectors with this structure. DFT-based precoders match in this case the channel structure providing a good beamforming gain. Note that the column vectors in a DFT codebook generate a grid of beams (GoB) that can be used to span all angular directions, as illustrated in Figure 5 for a planar array with 16 panels (subarrays), and $7\times 5$ antenna elements per subarray.

FIGURE 5

FIGURE 5. Illustration of the GoB generated with a 2D DFT codebook with an oversampling factor of 2. Orthogonal beams are marked in red.

In particular, 5G New Radio defines what is called a Type II DFT codebook, where an oversampling factor is introduced in the definition of array steering vectors for azimuth and elevation (Miao et al., 2018; Ericsson, 2017a, b; Samsung, 2017; Technical Specification Group Radio Access Network, 2017a, b). The goal of the oversampling factor, which is fixed and equal to 4 in 5G NR, is to mitigate the loss of beamforming gain in directions between orthogonal beams (straddling loss). This way, beam orthogonality is lost in general, and only beams spaced the oversampling factor are orthogonal. Figure 5 shows the set of orthogonal beams marked in red when an oversampling factor equal to 2 has been used in each dimension.

In our design, we will also consider a 2D oversampled DFT codebook for each subarray, although the oversampling factor O is not fixed, allowing the adjustment of the number of beams that illuminate footprints of different sizes for some given RF resources. The elements of the block diagonal RF precoding matrix ${\mathbf{F}}_{\text{RF}}\in {\mathrm{ℂ}}^{N^{\text{sat}}\times N_{\text{b}}}$ can be written in terms of this oversampled DFT codebook as

${\left[{\mathbf{F}}_{\text{RF}}\right]}_{n,I_{\text{RF}}\left(n\right)}=e^{-i\frac{2}{\lambda }\left({\left[{\mathbf{k}}_{\text{sat}}\right]}_{n,1}{\theta }_{I_{\text{RF}}\left(n\right),1}+{\left[{\mathbf{k}}_{\text{sat}}\right]}_{n,2}{\theta }_{I_{\text{RF}}\left(n\right),2}\right)},(14)$

for ${\left\{\left({\theta }_{m,1},{\theta }_{m,2}\right)\right\}}_m\subseteq {\left\{\frac{2\pi \left(n-1\right)}{O{N}_{\text{x}}^{\text{sub}}}\right\}}_n^{O{N}_{\text{x}}^{\text{sub}}}\times \{{\frac{2\pi \left(n-1\right)}{O{N}_{\text{y}}^{\text{sub}}}\}}_n^{O{N}_{\text{y}}^{\text{sub}}}$. Note that the m-th element of the codebook is contained in the non-zero entries of the mth column of ${\mathbf{F}}_{\text{RF}}$, that is ${\left[{\mathbf{F}}_{\text{RF}}\right]}_{n,m}=e^{-i\frac{2}{\lambda }\left({\left[{\mathbf{k}}_{\text{sat}}\right]}_{n,1}{\theta }_{m,1}+{\left[{\mathbf{k}}_{\text{sat}}\right]}_{n,2}{\theta }_{m,2}\right)}$. This expression is also proportional to the steering vector in the direction $\mathbf{v}\in {\mathcal{S}}_2$ with ${\left[\mathbf{v}\right]}_1=\frac{{\theta }_{m,1}}{\pi }$ and ${\left[\mathbf{v}\right]}_2=\frac{{\theta }_{m,2}}{\pi }$. To create the codebook, we prune the complete oversampled DFT dictionary, keeping only those beam-patterns pointing to the elliptical ROI. The oversampling factor can be used to adapt the RF resources to the size of the ROI and is the key parameter that drives the tesselation of the coverage area. By increasing O, the beam footprint is reduced, and the number of beampatterns pointing to the desired ROI increases. The optimal value for O, in terms of provided SNR, allows the creation of as many beams pointing to the ROI as available RF chains. Note that when increasing the oversampling factor, the dictionary provides a higher SNR over the ROI because the straddling loss decreases, but the interference between adjacent beams increases, leading to a reduction of the SINR. Therefore, the oversampling factor has to be carefully selected to find the appropriate tradeoff between SNR and SINR. The inter-beam interference problem can be addressed by designing a suitable digital precoding matrix ${\mathbf{F}}_{\text{BB}}$ that enables full frequency reuse, or by allocating different carriers to users served by adjacent beams, as proposed in (Angeletti and De Gaudenzi, 2020).

Once the codebook is created, the cells are defined as the areas inside the elliptical ROI where a given beam provides a higher gain than any other beam in the available codebooks for all the subarrays. Note that the DFT nature of the designed codebook creates a rectangular tiling of the ROI. Section 5 includes examples of the beam footprints when this type of codebook is assumed on board.

4 Coverage and Spectral Efficiency Calculation

4.1 Coverage

When operating with a transmitting beam-pattern $\mathbf{f}\in {\mathrm{ℂ}}^{N^{\text{sat}}}$ in the satellite beam codebook, and a receiving beam-pattern $\mathbf{w}\in {\mathrm{ℂ}}^{N^{\text{UT}}}$, the RSS is given by

$RSS\left[\text{W}\right]=|{\mathbf{w}}^{\text{H}}\mathbf{Hf}{|}^2=|\gamma {|}^2|<{\mathbf{a}}_{\text{u}}\left(\frac{{\mathbf{x}}_{\text{u}}-{\mathbf{x}}_{\text{sat}}}{\left|\left|{\mathbf{x}}_{\text{u}}-{\mathbf{x}}_{\text{sat}}\right|\right|}\right)+\frac{1}{\sqrt{K_{\text{r}}}}{\mathbf{a}}_{\text{R}},\mathbf{w}>{|}^2|<{\mathbf{a}}_{\text{sat}}^{\text{H}}\left(\frac{{\mathbf{x}}_{\text{sat}}-{\mathbf{x}}_{\text{u}}}{\Vert {\mathbf{x}}_{\text{sat}}-{\mathbf{x}}_{\text{u}}\Vert }\right),\mathrm{f}>{|}^2.(15)$

Note that this expression depends on the user and satellite positions in ${\mathbf{x}}_{\text{u}}$ and ${\mathbf{x}}_{\text{sat}}$ respectively. The different terms in Eq. 15 can be identified with the terms in the link budget equation. Thus $G_{\text{RX}}\left[\text{dB}\right]=10{\text{log}}_{10}\left(|<{\mathbf{a}}_{\text{u}}\left(\frac{{\mathbf{x}}_{\text{u}}-{\mathbf{x}}_{\text{sat}}}{\left|\left|{\mathbf{x}}_{\text{u}}-{\mathbf{x}}_{\text{sat}}\right|\right|}\right)+\frac{1}{\sqrt{K_{\text{r}}}}{\mathbf{a}}_{\text{R}},\mathbf{w}>{|}^2\right)$, the transmitter gain $G_{\text{TX}}\left[\text{dB}\right]=10{\text{log}}_{10}\left(|<{\mathbf{a}}_{\text{sat}}^{\text{H}}\left(\frac{{\mathbf{x}}_{\text{sat}}-{\mathbf{x}}_{\text{u}}}{\left|\left|{\mathbf{x}}_{\text{sat}}-{\mathbf{x}}_{\text{u}}\right|\right|}\right),\mathbf{f}>{|}^2\right)$, and $10{\text{log}}_{10}\left(|\gamma {|}^2\right)$ comprises all the losses $L{P}_{\text{cable}}\left[dB\right]+L{P}_{\text{at}}\left[dB\right]+L{P}_{\text{fs}}\left[dB\right]$.

To compute the combining gain at the UT side, we will consider the LoS direction to be known (in the system model we have assumed that the satellite position is known) and the Rician component of the channel to be unknown. Thus, the combiner ${\mathbf{w}}_{\text{u}}$ will be computed as the one maximizing $\frac{|\langle {\mathbf{a}}_{\text{u}}\left(\mathbf{v}\right),\mathbf{w}\rangle {|}^2}{{\left|\left|\mathbf{w}\right|\right|}^2}$ for the set of feasible $\mathrm{w}$ as discussed in Section 2.1. The average gain associated to this combiner can be computed as $\mathbb{E}\left(\frac{|\langle {\mathbf{a}}_{\text{u}}\left(\mathbf{v}\right)+\frac{1}{\sqrt{K_{\text{r}}}}{\mathbf{a}}_{\text{R}},{\mathbf{w}}_{\text{u}}\rangle {|}^2}{{\left|\left|{\mathbf{w}}_{\text{u}}\right|\right|}^2}\right)$. Since ${\mathrm{a}}_{\text{R}}$ is Gaussian and independent of any other variable, we can write the average gain as

$\mathbb{E}\left(\frac{{\left|\langle {\mathbf{a}}_{\text{u}}\left(\mathbf{v}\right)+\frac{1}{\sqrt{K_{\text{r}}}}{\mathbf{a}}_{\text{R}},{\mathbf{w}}_{\text{u}}\rangle \right|}^2}{{\left|\left|{\mathbf{w}}_{\text{u}}\right|\right|}^2}\right)=\mathbb{E}\left(\frac{{\left|\langle {\mathbf{a}}_{\text{u}}\left(\mathbf{v}\right),{\mathbf{w}}_{\text{u}}\rangle \right|}^2}{{\Vert {\mathbf{w}}_{\text{u}}\Vert }^2}\right)+\frac{1}{K_{\text{r}}}\mathbb{E}\left(\frac{{\left|\langle {\mathbf{a}}_{\text{R}},{\mathbf{w}}_{\text{u}}\rangle \right|}^2}{{\left|\left|{\mathbf{w}}_{\text{u}}\right|\right|}^2}\right).(16)$

Now, if we make use of the definition of ${\mathbf{w}}_{\text{u}}$ and consider the Rician component to be uniform, i.e. $\text{Σ}=\mathbf{I}$, we obtain the final expression for the receive gain as

$\mathbb{E}\left(\underset{\mathbf{w}}{\text{max}}\frac{{\left|\langle {\mathbf{a}}_{\text{u}}\left(\mathbf{v}\right),\mathbf{w}\rangle \right|}^2}{\Vert {\mathbf{w}}^2\Vert }\right)+\frac{1}{K_{\text{r}}}.(17)$

Since the LEO satellite is transmitting simultaneously through multiple beams with a gain $G_n^{\text{TX}}$, $n=1,\dots ,N_{\text{b}}$, besides considering the RSS or the SNR, we also need to understand the impact of interference when computing the coverage. To this aim, we will compute the interference as the RSS of the signals transmitted through undesired beam-patterns

$I\left[dW\right]=P_{\text{TX}}\left[dBW\right]-L{P}_{\text{cable}}\left[dB\right]+G_{\text{TX}}^{\text{interf}}\left[dB\right]-L{P}_{\text{at}}\left[dB\right]-L{P}_{\text{fs}}\left[dB\right]+G_{\text{RX}}\left[dB\right](18)$

with $G_{n_{\text{selected}}}^{\text{TX}}$ the gain corresponding to the desired beampattern and $G_{\text{TX}}^{\text{interf}}={\displaystyle \sum }_{n\ne n_{\text{selected}}}{G}_n^{\text{TX}}$. With these definitions, the signal to interference plus noise ratio (SINR) has the expression

$SINR=\frac{RSS\left[W\right]}{I\left[W\right]+{\sigma }^2\left[W\right]}.(19)$

4.2 Throughput and Spectral Efficiency

We consider the throughput expression given by the capacity $B_{\text{DL}}{\text{log}}_2\left(1+SINR\right)$ in bits per second. We also consider the throughput without interference $B_{\text{DL}}{\text{log}}_2\left(1+SNR\right)$ as a bound to what can be achieved with inter-beam interference mitigation techniques. Assuming that frequency division duplexing mechanisms enable beam sharing among close users illuminated by the same beam, we will also consider the spectral efficiency with interference ${\text{log}}_2\left(1+SINR\right)$ and without interference ${\text{log}}_2\left(1+SNR\right)$ as performance metric.

5 Simulation Results

Unless otherwise specified, we consider a LEO satellite orbiting at a 1,300 km height, covering an ellipsoidal ROI with semiradius 534.1 and 170.5 km. This area has been designed following the procedure described in 3.1 for a LEO constellation with 83 orbits with an inclination of ${53}^{\circ }$ and 53 satellites in each orbit. The downlink operates in the Ku band with a carrier frequency $f_{\text{DL}}=11.45\text{GHz}$ and a bandwidth $B_{\text{DL}}=250\text{MHz}$. A single polarization is assumed in the simulations.

To understand the role of the oversampling factor O in the beam codebook presented in Section 3.2, we show in Figure 6 different mesh designs covering the elliptical ROI for a $128\times 64$ satellite phased array and 8 RF chains in each dimension. Figure 6A shows the results when $O=2.49$, while Figure 6B shows the obtained rectangular mesh when $O=2$. Note that with $O=2$, only 39 out of 64 potential beams are used to illuminate the ROI, resulting in a very inefficient usage of the RF resources.

FIGURE 6

FIGURE 6. Examples of rectangular mesh covering an elliptical ROI when using the oversampled DFT beam codebook: (A) with $O=2.49$; (B) with $O=2$. Please note that this ROI is a toy example and does not correspond to the designed one for the simulations.

For the next simulations we will assume that the satellite phased arrays consist of $5\times 3$ sub-arrays of size $12\times 24$, with a total of $60\times 72$ antennas. Note that while we are arranging the sub-arrays into a larger plannar antenna array structure, each subarray will be transmitting an independent signal. We found out that, for this number of antennas and RF-chains, with an oversampling factor of 1.2 we can adjust all the available beams to the size of the ROI. The remaining system parameters are summarized in Table 1.

TABLE 1

TABLE 1. System parameters for the simulation of the downlink of a LEO SatCom system operating in the Ku band.

To understand the impact of different antennas at the user terminal, we will consider both the UPA and the leaky wave antenna (denoted as CRLH-LW) to have a number of elements equal to $24\times 24$. Beamforming with the planar antenna array at the user side introduces a gain of 27.6 dB. The gain provided by the leaky wave antenna is a function of the leakage factor, as illustrated in Figure 7A. The dependance of the antenna developed by Kymeta with the elevation angle (Kymeta, 2019) is also shown in Figure 7A.

FIGURE 7

FIGURE 7. Gain provided by the different antenna terminals: (A) Gain with perfect beam alignment. (B) Gain under misalignment. (C) Gain under phase quantization.

We consider now the effect of beam misalignment on the gain through Monte Carlo simulation for the phased antenna array and the leaky wave antenna. Figure 7B shows that, on average, the gains stay within reasonable levels, but the minimum values become critical under an error greater than $4\text{deg}$. The effects of phase shifter quantization for the phased array in a hybrid architecture using the beamforming method proposed in (Palacios et al., 2016) are shown in Figure 7C. It can be seen that a hybrid architecture offers a much better performance than an analog one under phase quantization, but both perform very similarly when the phase shifters have a high resolution.

Next, we analyze the coverage of the elliptical ROI provided by the beam codebook at the satellite when the UPA is considered as antenna terminal. As shown in Section 4.1, $L{P}_{\text{fs}}$, $G_{\text{TX}}$ and $G_{\text{TX}}^{\text{interf}}$ depend on the user location. The coverage provided by the satellite’s codebook in terms of SNR and SINR can be seen in Figure 8. While the SNR shows a good coverage of the ROI for an oversampling factor of 1.2, when visualizing the SINR we can observe areas with high interference. This is a consequence of using the DFT-type codebook employed in 5G NR, since its associated rectangular mesh creates points that are covered with the same gain by 4 different satellite beams. When selecting one of the beams, the other three will create a strong interference, leading to a SINR not larger than $-4.77\left[\text{dB}\right]$. Furthermore, the design of static codebooks will always lead to points in which at least 3 beam-patterns have equal gain, thus leading to a SINR never larger than $-3\left[\text{dB}\right]$.

FIGURE 8

FIGURE 8. (A) Cells associated to the different beampatterns; point marked in green show the location where the beampattern gain is maximum. (B) SNR when using the best beam for each cell. (C) SINR when using the best beam for each cell.

This SNR and SINR translate into the spectral efficiency shown in Figure 9. Note that there is a significant inter-beam interference when operating with DFT-type codebooks as proposed in LTE and 5G NR. Future work needs to be developed to design alternative analog codebooks and hybrid beamforming strategies at the satellite side to reduce inter-beam interference.

FIGURE 9

FIGURE 9. (A) Spectral efficiency when ignoring interference. (B) Spectral efficiency when considering interference.

To further extend the coverage analysis, in Figure 10 we see how the SNR, SINR and their corresponding spectral efficiencies are distributed over the ROI as CDF curves. It is easy to appreciate the huge impact of inter-beam interference, hinting how important the design of good interference mitigation techniques is for satellite communications.

FIGURE 10

FIGURE 10. Cumulative distribution function inside the ROI of: (A) SNR and SINR. (B) Spectral efficiency computed for the SNR and SINR.

Figure 11 shows the satellite’s angular speed relative to the user and the Doppler effect due to the satellite’s movement as a function of the user’s relative location. The relative angular speed shows significantly small values. This means that the beampattern switching requirements of a given antenna design that enable tracking of the satellite’s position for beamforming are not particularly challenging. Note that this holds when tracking the satellite direction, not its corresponding azimuth and elevation angles. If the angular relative speed is translated into a maximum azimuth and elevation angular variation with time, the azimuth angle can present sharp changes, as illustrated in Section. 2.2.2.

FIGURE 11

FIGURE 11. (A–C) Doppler effect for the central frequency caused by the satellite’s movement, Earth rotation and UT movement respectively. (D) Satellite’s angular relative speed with respect to the user terminal.

Finally, we analyze the SNR degradation over time due to the satellite’s movement for a user initially located at the center of the elliptical ROI. We have considered a beam selection strategy based on maximizing SNR. Figure 12A shows an example of different satellite positions and associated selected beams. Figure 12B shows the SNR corresponding to the different satellite positions illustrated in (a). It can be concluded that, for this codebook design and system parameters, the user is illuminated by the same beam for a period of around 30 s. After that, another beam has to be selected to keep a proper gain. Note that even after updating the best beam, there is a non-negligible gain loss due to the satellite movement, as also illustrated in Figure 12B.

FIGURE 12

FIGURE 12. (A) Example of different satellite positions and associated beams. (B) SNR degradation over time due to the satellite’s movement.

6 Conclusion

Research on integrating LEO-based satellite networks with terrestrial networks is still in its infancy. We implemented a complete simulator of the downlink of a massive MIMO LEO SatCom system to understand the physical layer challenges and evaluate a hybrid beamforming strategy based on an analog codebook. We designed the satellite footprint as an elliptical ROI given the number of orbital planes and number of satellites per orbit. We also designed a static DFT-type codebook, as commonly used in current cellular standards, to adjust the beamwidth and RF resources to the size of the satellite footprint by means of an interpolation factor. Simulation results show the relevance of the antenna terminal design to achieve an acceptable SNR. Our numerical experiments also show that DFT codebooks introduce a high inter-beam interference which reduces the ability of the system to provide a high spectral efficiency at all locations inside the ROI. Designing fully or partially connected hybrid beamforming architectures with per antenna power constraints and a digital precoder that mitigates inter-beam interference is an open research challenge that needs to be addressed to overcome the SINR performance of current approaches. The design of low complexity dynamic codebooks that better adapt the resources to the specific location of the users and reduce the SNR loss due to the satellite movement is a critical challenge. Our simulator also shows that the variation of the relative angular speed between the satellite and the UT is smooth when working with spherical coordinates. This leads to the need of designing beam tracking methods for the UT that track the satellite position rather than its variation in elevation and azimuth. Our numerical experiments also show the effectiveness of exploiting satellite position information at the UT side to avoid the need of channel estimation and tracking, which will provide old estimates due to the long transmission delay. Finally, an accurate model for the error in the satellite’s position and trajectory is also necessary to better understand its impact in the performance of the satellite beam tracking strategies.

Data Availability Statement

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

Author Contributions

JP and NG-P contributed to the design and development of the investigation, including mathematical derivations and code for the simulations. NG-P and CM provided technical direction. JP and NG-P wrote the draft of the manuscript. CM, TS and C-HW participated in technical discussions of the results and provided comments of the manuscript.

Funding

This work has been partially funded by Toyota Motor North America.

Conflict of Interest

Authors TS and C-HW were employed by the company Toyota Motor North America.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Angeletti, P., and De Gaudenzi, R. (2020). A Pragmatic Approach to Massive MIMO for Broadband Communication Satellites. IEEE Access 8, 132212–132236. doi:10.1109/access.2020.3009850

CrossRef Full Text | Google Scholar

De Gaudenzi, R., Angeletti, P., Petrolati, D., and Re, E. (2020). Future Technologies for Very High Throughput Satellite Systems. Int. J. Satell Commun. Netw. 38, 141–161. doi:10.1002/sat.1327

CrossRef Full Text | Google Scholar

del Portillo, I., Cameron, B. G., and Crawley, E. F. (2019). A Technical Comparison of Three Low Earth Orbit Satellite Constellation Systems to Provide Global Broadband. Acta Astronautica 159, 123–135. doi:10.1016/j.actaastro.2019.03.040

CrossRef Full Text | Google Scholar

Devillers, B., and Pérez-Neira, A. (2011). “Advanced Interference Mitigation Techniques for the Forward Link of Multi-Beam Broadband Satellite Systems,” in 2011 Conference Record of the Forty Fifth Asilomar Conference on Signals, Systems and Computers (ASILOMAR), 1810–1814. doi:10.1109/ACSSC.2011.6190334

CrossRef Full Text | Google Scholar

[Dataset] Ericsson, (2017a). On Basis Design for Type I and Type II Codebooks.

Ericsson, (2017b). Type II CSI for Beamformed CSI-RS and Hybrid Operation. Tech. Rep. R1-1708699. Hangzhou, China. Meeting RAN1#89, 3GPP.

Kim, J., Yun, M. Y., You, D., and Lee, M. (2020). “Beam Management for 5G Satellite Systems Based on NR,” in 2020 International Conference on Information Networking (ICOIN), 32–34. doi:10.1109/icoin48656.2020.9016571

CrossRef Full Text | Google Scholar

Kodheli, O., Guidotti, A., and Vanelli-Coralli, A. (2017). “Integration of Satellites in 5G through LEO Constellations,” in GLOBECOM 2017 - 2017 IEEE Global Communications Conference, 1–6. doi:10.1109/glocom.2017.8255103

CrossRef Full Text | Google Scholar

[Dataset] Kymeta (2019). Link Budget Calculations for A Satellite Link with an Electronically Steerable Antenna Terminal.

Mehdipour, A., and Eleftheriades, G. V. (2014). Leaky-Wave Antennas Using Negative-Refractive-Index Transmission-Line Metamaterial Supercells. IEEE Trans. Antennas Propagat. 62, 3929–3942. doi:10.1109/TAP.2014.2322882

CrossRef Full Text | Google Scholar

Mehdipour, A., Sazegar, M., and Stevenson, R. (2019). “Broadband WAIM Metasurface Structure for Electronically Beam Scanning Holographic Antenna for Ku-Band Satellite Communications,” in 2019 IEEE International Symposium on Antennas and Propagation and USNC-URSI Radio Science Meeting, 429–430. doi:10.1109/APUSNCURSINRSM.2019.8888376

CrossRef Full Text | Google Scholar

Miao, H., Mueck, M. D., and Faerber, M. (2018). “Amplitude Quantization for Type-2 Codebook Based CSI Feedback in New Radio System,” in 2018 European Conference on Networks and Communications (EuCNC). 1–9. doi:10.1109/EuCNC.2018.8442609

CrossRef Full Text | Google Scholar

Palacios, J., De Donno, D., Giustiniano, D., and Widmer, J. (2016). “Speeding up mmWave Beam Training through Low-Complexity Hybrid Transceivers,” in 2016 IEEE 27th Annual International Symposium on Personal, Indoor, and Mobile Radio Communications (PIMRC), 1–7. doi:10.1109/PIMRC.2016.7794709

CrossRef Full Text | Google Scholar

Samsung (2017). WF on Type I and II CSI Codebooks. Tech. Rep. R1-1709232. Hangzhou, China. Meeting RAN1#89, 3GPP.

Sikri, D., and Jayasuriya, R. M. (2019). Multi-Beam Phased Array with Full Digital Beamforming for SATCOM and 5G. Microwave J. 60 (4), 64–79.

Google Scholar

[Dataset] Technical Specification Group Radio Access Network (2017a). TR38.211. NR; Physical Channels and Modulation (Release 15).

[Dataset] Technical Specification Group Radio Access Network (2017b). TR38.214. NRPhysical Layer Procedure for Data (Release 15).

[Dataset] Technical Specification Group Radio Access Network (2019). TR38.821. Solutions for NR to Support Non-terrestrial Networks (NTN) (Release 16).

Xia, S., Jiang, Q., Zou, C., and Li, G. (2019). Beam Coverage Comparison of LEO Satellite Systems Based on User Diversification. IEEE Access 7, 181656–181667. doi:10.1109/access.2019.2959824

CrossRef Full Text | Google Scholar

You, L., Li, K.-X., Wang, J., Gao, X., Xia, X.-G., and Ottersten, B. (2020). Massive MIMO Transmission for LEO Satellite Communications. IEEE J. Select. Areas Commun. 38, 1851–1865. doi:10.1109/jsac.2020.3000803

CrossRef Full Text | Google Scholar