Analyzing Downlink Coverage in Clustered Low Earth Orbit Satellite Constellations: A Stochastic Geometry Approach

Satellite communications have garnered interest as a promising solution to achieve ubiquitous connectivity. According to the International Telecommunication Union report [1], nearly half of the global population remains unconnected. Satellites are expected to address this connectivity gap, offering a viable solution that alleviates the cost burden faced by telecom operators when deploying terrestrial base stations [2]. Traditionally, geostationary orbit (GEO) satellites were the primary choice for satellite networks due to the extensive coverage facilitated by their high altitudes. However, the extremely high altitudes of GEO satellites present challenges such as high latency and low area spectral efficiency, making them less suitable for 5G applications. The industry has shifted towards adopting low Earth orbit (LEO) satellites to reduce latency, increase service density, and facilitate cost-effective launches [3, 4, 5]. Companies like SpaceX, OneWeb, and Amazon are planning to deploy large constellations in LEOs, often referred to as mega-constellations [6]. However, in mega-constellations, launching additional satellites does not necessarily increase throughput and could even lead to performance degradation due to severe inter-satellite interference. To tackle this challenge, satellite cooperation concepts such as LEO remote sensing [7] and satellite clusters [8],[9] have been introduced to manage satellite networks effectively and mitigate inter-satellite interference.

As mega-constellations complicate simulation-based network performance analyses, the need for new tools arises for effective assessment. Stochastic geometry is introduced to mathematically analyze performance, modeling the random behavior of nodes in wireless networks, including base stations and users, through point processes such as binomial point processes (BPPs) or Poisson point processes (PPPs) [10]. Previous studies [11, 12, 13, 14, 15, 16] have conducted stochastic geometry-based analyses for wireless networks using various performance metrics for terrestrial networks. A framework for modeling wireless networks was proposed in [11], addressing the outage probability, network throughput, and capacity. In [12], the distribution of base stations modeled by a PPP was compared to the actual distribution in cellular networks, and the coverage probability and average rate were evaluated. The authors in [13] modeled downlink heterogeneous cellular networks with the PPP and analyzed the coverage probability and average rate. In [14], the coverage probability for uplink cellular networks modeled by the PPP was studied. The authors in [15] derived the downlink coverage probability considering base station cooperation within the PPP framework. A BPP model for a cache-enabled networks was used to characterize the downlink coverage probability and network spectral efficiency in [16].

Building on stochastic geometry-based analyses in terrestrial networks, recent works have extended its application to satellite networks. Given the envisaged prevalence of large-scale satellite constellations, it is crucial to analyze the network performance in scenarios where multiple satellites serve a user, moving beyond a single-satellite setup. In this work, we focus on satellite cluster networks, a pertinent solution in mega-constellations, and employ stochastic geometry to characterize the coverage probability, providing an effective means of estimating the impact of system parameters.

We consider a downlink satellite network depicted in Fig. 1. Satellites are located on a sphere with radius $R_{\mathrm{S}}$ according to a homogeneous spherical Poisson point process (SPPP) $\Pi=\{\mathbf{x}_{1},\ldots,\mathbf{x}_{\mathsf{N}_{\mathrm{S}}}\}$ where $\mathbf{x}_{l}$, $l\in\{1,\ldots,\mathsf{N}_{\mathrm{S}}\}$, is the position of the satellite $l$ in a spherical coordinate system with elevation and azimuth angles spanning from $0$ to $\pi$ and $0$ to $2\pi$, respectively. The number of total satellites, ${\mathsf{N}_{\mathrm{S}}}$, follows the Poisson distribution with density $\lambda_{\mathrm{S}}$. Single-antenna users, positioned on the Earth’s surface with radius $R_{\mathrm{E}}$, are also distributed according to a homogeneous SPPP, denoted as $\Pi_{\mathrm{U}}=\{\mathbf{u}_{1},\ldots,\mathbf{u}_{\mathsf{N}_{\mathrm{U}}}\}$, where ${\mathsf{N}}_{\mathrm{U}}$ is the number of total users on the Earth, which follows the Poisson distribution with density $\lambda_{\mathrm{U}}$. According to the Slivnyak’s theorem [10], we focus on a typical user located at $\mathbf{u}_{1}=(0,0,R_{\mathrm{E}})$ for our downlink performance analyses without loss of generality.

To derive the coverage probability, it is necessary to establish a basic performance metric and mathematical preliminaries. We assume an interference-limited network and employ the signal-to-interference ratio (SIR) as a performance metric. This choice encapsulates the scenario where interference becomes predominant with the increasing number of satellites, making noise negligible. We also approximate the accumulated received power from the satellites within and outside the cluster area as Gamma random variables, consolidating these compound terms into their respective single random variables.

In this section, we propose the coverage probability of the satellite cluster network, with a focus on downlink satellite cooperation in the cluster area. We utilize the SIR that incorporates joint transmissions in the cluster area, as defined in Section III-A, and apply the Gamma approximation outlined in Section III-B to derive the lower and upper bounds for the coverage probability. We assume the mutual independence of the processes governing the distribution of satellites within and outside the cluster area, given their distinct spatial regions.

The coverage probability in the interference-limited regime is the probability that the SIR at the typical user is greater than or equal to a threshold $\gamma$. This probability is derived as

where $\gamma$ is the SIR threshold. To derive the coverage probability based on (14), we use the parameters given in (10) and (11) to consider the CDF of the approximated Gamma random variable $\tilde{\mathsf{I}}$. Similarly, when obtaining the coverage probability based on (15), we utilize the parameters provided in (12) and (13) to incorporate the complementary CDF of the corresponding approximated Gamma random variable $\tilde{\mathsf{D}}$. For the remainder of this section, we opt for (14) to attain the coverage probability with the approximated Gamma random variable $\tilde{\mathsf{I}}$.

In order to obtain the coverage probability as per (14), it is necessary to compute the Laplace transform of $\mathsf{D}$. This process commences with obtaining the conditional Laplace transform of $\mathsf{D}$, which is then extended to the unconditional case through marginalization over the Poisson distribution. Before obtaining the conditional Laplace transform, we introduce the PDF of the distance between the typical user and a satellite located in the observable spherical dome $\mathcal{A}$ in the following lemma.

From the result, we derive the conditional Laplace transform of $\mathsf{D}$. Let $\mathsf{N}_{\mathrm{S,clu}}$ denote the number of satellites in the cluster area $\mathcal{A}_{\mathrm{clu}}$. Fixing the number of cooperating satellites at $\mathsf{N}_{\mathrm{S,clu}}=L$, the locations of the $L$ satellites can be modeled using a BPP [34]. Then the conditional Laplace transform is given in the following lemma.

In solving the integral in the conditional Laplace transform, computation complexity may arise. However, assuming the Rayleigh fading for satellite channels, i.e., $m=1$, reduces the integral to a closed-form solution. The simplified conditional Laplace transform for $m=1$ is obtained in the following corollary.

By marginalizing the conditional Laplace transform over $L$, the Laplace transform of the accumulated cluster power $\mathsf{D}$ is obtained as

where (a) is obtained using the Poisson distribution, and (b) is determined through the Taylor series expansion of the exponential function, $e^{x}=\sum_{n=0}^{\infty}\frac{x^{n}}{n!}$. It is noteworthy that the conditional Laplace transform in (IV) becomes independent of $L$ through the $\frac{1}{L}$th root operation on (IV).

Due to the non-integer property of the shape parameters $k_{\tilde{\mathsf{I}}}$ and $k_{\tilde{\mathsf{D}}}$ in most cases, obtaining the exact form of the coverage probability entails the utilization of the Gamma random variable’s CDF, which involves infinite summations. These infinite operations impose computational complexity and render the mathematical analysis of the results challenging. To release this issue without compromising generality, we proceed to establish the lower and upper bounds of the coverage probability using the Erlang distribution while considering a finite number of terms $\tilde{k}$ in both differentiation and summation [36]. Using the definition of the coverage probability (14) and the Laplace transform derived (IV), we establish the coverage probability bounds in the following theorem. It is noteworthy that the lower and upper bounds, which will be derived shortly, coincide when the shape parameter, $k_{\tilde{\mathsf{I}}}$ or $k_{\tilde{\mathsf{D}}}$, is an integer.

The coverage probability based on (15) then incorporates the approximated Gamma random variable $\tilde{\mathsf{D}}$ and the Laplace transform of $\mathsf{I}$. More detailed results are presented in the following theorem.

From the above results, we can see that the coverage probability depends on various system parameters, including threshold $\gamma$, satellite density $\lambda_{\mathrm{S}}$, satellite altitude $R_{\mathrm{S}}$, cluster area determined by $\phi_{\mathrm{clu}}$, path loss exponent $\alpha$, and Nakagami parameter $m$. The primary distinction between adopting $\tilde{\mathsf{I}}$ or $\tilde{\mathsf{D}}$ pertains to their respective shape parameters, which are determined by the difference in distances, as shown in (10) and (12). Notably, the shape parameter for $\tilde{\mathsf{I}}$ is usually higher than that of $\tilde{\mathsf{D}}$ because the difference between $R_{\mathrm{clu}}$ and $R_{\mathrm{max}}$ is much larger than the difference between $R_{\mathrm{min}}$ and $R_{\mathrm{clu}}$ particularly when considering a small cluster polar angle $\phi_{\mathrm{clu}}$. As shown, $\tilde{k}$ determines the extent of differentiation and summation required in the coverage probability calculation. The large shape parameter of $\tilde{\mathsf{I}}$ then narrows the gap between the lower and upper bounds of the coverage probability but increases the computational burden. Conversely, employing a small shape parameter with $\tilde{\mathsf{D}}$ can potentially result in a wider gap between the two bounds compared to $\tilde{\mathsf{I}}$, but it enables the analyses of network with a massive number of satellites while reducing the computational load. More details about this trade-off can be found in simulation results.

We further simplify the derivative of the Laplace transform $\mathcal{L}_{\{\mathsf{D}\}}(s)$ leveraging the differentiation property of the exponential function. Utilizing Faà di Bruno’s formula [37], the $n$-th derivative of $\mathcal{L}_{\{\mathsf{D}\}}(s)$ is expressed as

where the Bell polynomial is defined as ${B_{n,q}(x_{1},\ldots,x_{n-q+1})}$ ${=\sum\frac{n!}{a_{1}!\cdots a_{n-q+1}!}\left(\frac{x_{1}}{1!}\right)^{a_{1}}\cdots\left(\frac{x_{n-q+1}}{(n-q+1)!}\right)^{a_{n-q+1}}}$, and the sum is over all sequences $a_{1},a_{2},\ldots,a_{n-q+1}$ of non-negative integers such that ${a_{1}+a_{2}+\cdots+a_{n-q+1}=q}$. Then, we can rewrite the lower and upper bounds of the coverage probability as follows

As can be seen in (IV), differentiating the exponent of the Laplace transform is now required to obtain the bounds of the coverage probability. For the special case under the Rayleigh fading, i.e., $m=1$, we simplify the $n$-th derivative of the exponent of $\mathcal{L}^{m=1}_{\{\mathsf{D}\}}(s)$ in the following corollary.

In this section, we validate the Gamma approximation and the derived coverage probability bounds through Monte Carlo simulations. We also evaluate the impact of the cluster area on the coverage probability.

In our numerical experiments, we set the satellite altitude to $R_{\mathrm{S}}-R_{\mathrm{E}}=500$ km and restrict satellite visibility to $\theta_{\mathrm{min}}=25^{\circ}$, with a defined polar angle of the cluster area $\phi_{\mathrm{clu}}=1.6^{\circ}$. Utilizing a sectored antenna gain, we employ the normalized antenna gain $G_{\mathrm{o}}^{\mathrm{t}}/G_{\mathrm{i}}^{\mathrm{t}}=0.1$, indicating a 10 dB gain advantage for satellites in the cluster area compared to those outside. For the coverage probability calculations, we explore two scenarios: one with an average of $\lambda_{\mathrm{S}}|\mathcal{A}|=50$ satellites in the cluster area and the other with $\lambda_{\mathrm{S}}|\mathcal{A}|=300$. With these parameters, the average number of satellites on the satellite sphere, whose surface area is calculated as $4\pi R_{\mathrm{S}}^{2}$, is at least 10k, and that on the cluster area exceeds two. Specific satellite numbers for the scenarios are detailed in Table I. In addition, the first scenario employs the distribution of $\tilde{\mathsf{I}}$, the approximated Gamma random variable for the accumulated interference power $\mathsf{I}$, while the second scenario utilizes the distribution of $\tilde{\mathsf{D}}$, the approximated Gamma random variable for the accumulated cluster power $\mathsf{D}$. The upcoming results will offer a more detailed explanation, demonstrating how the utilization of these random variables enables achieving either elaborateness or computational efficiency.

In this paper, we proposed the mathematical analyses of the performance of satellite cluster networks, where satellites in the cluster area collaborate to serve users, particularly in mega-constellations. To establish the lower and upper bounds of the coverage probability, we first derived the key parameters of the approximated Gamma random variables to handle the compound terms in the SIR and compared the advantages of each in terms of tightness and complexity. Leveraging the distribution of the approximated Gamma random variable, we suggested the lower and upper bounds of the coverage probability, which depend on system parameters. While the results in this work relied on random variable approximations and provided lower and upper bounds, these bounds effectively showed the network’s performance with sufficient tightness to simulation results. Moreover, our analyses not only reduced computational burden through finite operations but also expected to offer insights into the impact of system parameters on the coverage probability in satellite cluster networks.