## I Introduction

The emerging Internet-of-Things (IoT) enables information exchange among physical objects, and enlightens the future development for a wide range of applications [1, 2, 3], such as smart metering, e-health, fleet management, smart cities, etc. Therefore, supporting a variety of IoT applications becomes one major task for future wireless communication systems. Considering the engineering and operation cost, existing cellular communication networks are mainly deployed at populated areas, with the aim to meet the need for human-to-human (H2H) communications. However, compared with H2H devices, IoT devices are widely deployed in remote areas like deserts, coastal waters, and forests [4]. Therefore, it is challenging and cost-inefficient to provide services for remote IoT devices with existing cellular communication infrastructures.

As an alternative to cellular communication networks, the satellite communication network provides a promising solution to supporting IoT services [5]. Specifically, according to recent reports, thousands of low-earth orbit (LEO) satellites are planned to be launched by SpaceX [7] and OneWeb [8]. Driven by these advances, a dense constellation of LEO satellites will be established to provide seamless Internet coverage for terrestrial users. Compared with existing cellular communication network and traditional geostationary earth orbit (GEO) satellite network, the LEO satellite network exhibits the following advantages in enabling IoT applications [6]:

(i) Compared with existing cellular network, the terrestrial-satellite link (TSL) in LEO satellite communication is more robust in different terrestrial environment. Therefore, the LEO satellite network is more tolerant to geological disasters and extreme topographies like cliffs, valleys, and steep slopes.

(ii) Compared with cellular network, the LEO satellite network requires less support from terrestrial infrastructures such as base stations. Therefore, for IoT devices deployed in remote areas like deserts, forests and oceans, the LEO satellite network is more cost-efficient. More importantly, the LEO satellite network is also preferred for achieving global coverage due to its dense constellation.

(iii) The orbital altitude of LEO satellites (normally several hundred kilometers) is much lower than that of GEO satellites (approximately ). Therefore, the propagation delay and path loss in LEO satellite communications are smaller than those in GEO satellite communications [4]. Furthermore, lower transmission power is required to reach LEO satellites than to reach GEO satellites, which is preferable for IoT devices with stringent requirements on power consumption.

(iv) In GEO satellite communications, terrestrial devices can only access the satellite with fixed elevation angle. Therefore, GEO satellite communications is vulnerable to obstacles between the devices and the satellite receiver. However, in LEO satellite communications, IoT devices can access the mobile LEO satellites with flexible elevation angles. Therefore, LEO satellite-enabled IoT is tolerant to terrestrial obstacles.

### I-a Related Works

Motivated by above-mentioned advantages, some research has been conducted to exploit LEO satellites for future wireless communication network. For example, a dense LEO satellite access network (LEO-SAN) was proposed in [9], which integrates terrestrial communications with satellite communications. In this LEO-SAN, different physical-layer techniques such as interference management, diversity techniques, and cognitive radio schemes were investigated to achieve seamless coverage. A similar network architecture was proposed in [10] to integrate terrestrial-satellite network (TSN) into 5G and beyond to achieve efficient data offloading. The application of power-domain non-orthogonal multiple access (NOMA) was investigated in [11] for various satellite architectures. For IoT applications, an overview was presented in [4] on satellite-enabled Internet of Remote Things (IoRT), where various enabling techniques are discussed, such as the MAC protocol, resource allocation and transmission management. Furthermore, other topics including the constellation structure, spectrum allocation, and routing protocols are discussed in [6]. The inter-plane inter-satellite link (ISL) was considered in [12] to further extend the coverage of IoT applications.

According to the service type [13, 14]

, IoT devices are intermittently activated with a certain probability and short data packets, i.e., active IoT devices perform random access (RA) for sporadic data transmission. Conventional RA schemes for satellite communications are mainly based on slotted ALOHA protocols

[4, 15, 16]. In addition, a divide-and-conquer scheme was proposed in [17] to allocate time slots to terminals based on service demands. Generally, these RA schemes assume sufficient access resources and static TSL, which becomes infeasible for IoT applications with massive connectivity and LEO satellites with high mobility.Confronted with the massive connectivity in IoT and the rapidly-changing TSL to mobile LEO satellites, grant-free RA (GF-RA) schemes are preferred due to its spectral efficiency and low access delay. In GF-RA schemes, activated devices share the same access resource and directly transmit their data packets (along with pilot sequences), without applying for the grant from the satellite receiver. In this way, the signaling overhead can be reduced, which improves the transmission efficiency for the short data packets of IoT devices.

Inspired by these advantages, a space diversity-based GF-RA scheme was proposed for satellite-enabled IoT [18]. However, one crucial task for the satellite receiver, i.e. the joint user activity detection (UAD) and channel estimation (CE) was bypassed in [18] by roughly modeling the TSL as an erasure-collision channel. To solve this joint UAD and CE problem, a typical solution is to formulate this task as a compressed sensing (CS) problem. In [19, 20, 21, 22]

, the pilot sequences from different devices serve as the sensing matrix in this CS problem, and the approximate message passing (AMP) algorithm was proposed for different system models. In addition, the variational Bayesian inference was employed in

[23, 24, 25], where the UAD and CE are solved by the mean-field message passing algorithm [26] and the Gaussian message passing algorithm [27, 28, 29], respectively.### I-B Motivations

Existing RA schemes for satellite communications [15, 16, 17] rely on static TSL and sufficient access resources. However, in LEO satellite-enabled IoT, the TSL is rapidly changing due to the mobility of the LEO satellite, and the massive connectivity from IoT devices may cause shortage of access resources. Therefore, existing RA schemes for satellite communications fail to match with the LEO satellite-enabled IoT. Although GF-RA schemes can improve the spectral efficiency and lower the transmission delay, the joint UAD and CE algorithms in existing GF-RA schemes [19, 20, 21, 22, 23, 24, 25] are designed for the Rayleigh channel, instead of the TSL channel. In the Rayleigh channel model, the line-of-sight (LoS) component in wave propagation is negligible, while the scattering component is dominant. However, both the LoS component and the scattering components should be considered for the TSL. To the best of our knowledge, there is little work tailored for the GF-RA scheme in LEO satellite-enabled IoT, which motivates our research in this paper.

### I-C Contributions

In this paper, we propose a terrestrial-satellite GF-RA scheme for the LEO satellite-enabled IoT. In order to address the joint UAD and CE problem, we propose a Bernoulli-Rician message passing with expectation maximization (BR-MP-EM) algorithm. The iterative message passing process of this BR-MP-EM algorithm is described on a factor graph and can be divided into two stages, i.e., the inner iterations and outer iterations. In the inner iterations, the Bernoulli messages and Rician messages are jointly updated for the UAD and CE. The channel impairments, i.e. the random propagation fading and phase shift, are treated as hyper-parameters. Then the EM method is employed in the outer iterations to update these hyper-parameters. The major contributions of this paper are listed as follows:

(i) A terrestrial-satellite GF-RA scheme is proposed for LEO satellite-enabled IoT, where active devices share the same access resource and directly transmit their data packets. This GF-RA scheme could avoid excessive signaling overhead, and thus improve the transmission efficiency for the short data packets of IoT devices.

(ii) With a typical channel model for the TSL, we derive the joint Bernoulli-Rician message passing in the inner iterations of the BR-MP-EM algorithm. The Bernoulli message update is derived for the UAD, while the Rician message update is derived for the CE.

(iii) In the outer iterations of the BR-MP-EM algorithm, the EM update is derived to estimate the channel impairment-related hyper-parameters, based on the output of the inner iterations. This EM update is independent from specific distributions of the hyper-parameters, and therefore robust against unknown channel impairments.

### I-D Potential Applications

A terrestrial-satellite GF-RA scheme is proposed in this paper, and we derive the BR-MP-EM algorithm to solve the joint UAD and CE problem in this GF-RA scheme. The proposed GF-RA scheme enables sporadic transmission from IoT devices to LEO satellites, and can be potentially applied to a wide range of IoT applications, especially when the IoT devices are remotely deployed so that the assistance of LEO satellites becomes essential. Some typical examples are smart grid, emergency management, inshore and offshore windmills, land environment monitoring and ocean monitoring [4].

### I-E Section Organization

This paper is organized as follows. In Section II, the system model is explained, and the joint UAD and CE problem is formulated for the terrestrial-satellite GF-RA system. In Section III, the BR-MP-EM algorithm is proposed, with detailed derivations on the Bernoulli-Rician message passing and the EM update. Simulation results are presented in Section IV to evaluate the UAD and CE performance of the proposed BR-MP-EM algorithm, as well as its robustness against unknown channel impairments. Finally, Section V concludes this paper.

## Ii System Model and Problem Formulation

### Ii-a System Model

As illustrated in Fig. 1, we consider a terrestrial-satellite GF-RA system in LEO satellite-enabled IoT. Within the coverage area of a serving satellite, there are potential devices, each of which is randomly activated with a probability . In each round of GF-RA, all the activated devices share the same access resource and transmit their data packets to the serving satellite (along with their unique pilot sequences). Considering the characteristics of dense LEO constellation and IoT applications, we adopt the following assumptions for this terrestrial-satellite GF-RA system:

(i) Each device in this system refers to one physical access point (AP), which is equipped with a dedicated terrestrial-satellite terminal (TST) to facilitate transmission to the LEO satellite [9, 10]. We assume that each AP provides coverage for a small number of nearby IoT terminals over wired or wireless terrestrial link. In this way, power-constrained sensors and IoT terminals can access the LEO satellite via the AP. Due to the low activity of these IoT terminals, the AP is also randomly activated in each round of GF-RA with a probability.

(ii) We assume that all these devices are stationary, i.e., the devices are deployed at fixed locations with no mobility [4]. Furthermore, the surrounding terrestrial environment of these devices is assumed sufficiently static.

(iii) Since LEO satellites orbit the earth with fixed time period [6], we assume that aided by calibration from terrestrial control stations, the satellite receiver can maintain time synchronization with terrestrial devices. In addition, the timing advance information for each device’s transmission is configured according to its location. In this way, data packets from different devices are aligned in time at the receiver.

(iv) On top of assumptions (ii) and (iii), we assume that activated devices perform GF-RA when the serving satellite appears at a fixed orbital position. In this way, in different rounds of GF-RA, each device can perform transmission at approximately fixed elevation angle.

For the TSL in this terrestrial-satellite GF-RA system, we adopt a widely-applied channel model, which is shown to match well with measured land mobile satellite (LMS) channel data [30]. Specifically, the channel from each device to the LEO satellite is modeled as [9, 11, 30, 31, 32, 33]

(1) |

where the combined effect of the scattering components between each device and the serving satellite (denoted by dashed lines in Fig. 1) is characterized by a Rayleigh-distributed amplitude and a uniformly random phase . On the other hand, the LoS component between each device and the serving satellite (denoted by solid lines in Fig. 1) is characterized by a Nakagami-distributed amplitude and a deterministic phase . It is noted that the random amplitude is mainly caused by the shadowing effect [30].

On top of the LMS channel model (1), we add two additional factors to account for the channel impairments, i.e., the random propagation fading and phase shift , which will be explained later. Then, the terrestrial-satellite channel for each device is modeled as

(2) |

Compared with (1), a channel impairment term is included in (2), and we use subscript in the notations to distinguish the channel and wave-propagation components for different devices. It is noted that we can describe the scattering component in (2

) by a complex Gaussian distribution with variance

, i.e., . According to assumptions (ii) and (iv), the stationary devices perform GF-RA at fixed elevation angles. As a result, the shadowing effect on the LoS component is also sufficiently static in each round of GF-RA, so that the LoS amplitude can be considered as a constant in the proposed GF-RA scheme. In addition, we assume that all the devices experience the same phase shift and propagation fading in each round of GF-RA, and the reason for this assumption is explained in Remark 1.### Ii-B Problem Formulation

Assume that each device has a unique pilot sequence with length . When device is activated, is transmitted along with its data packet to facilitate the joint UAD and CE at the satellite receiver. Then, the -th received pilot symbol is

(3) |

where is the set of devices, is the activity indicator, i.e., if device is activated. Otherwise, . is the terrestrial-satellite channel of device in (2), and is the additive white Gaussian noise (AWGN) with variance

. Considering all the received pilot symbols, we have the following received pilot vector

(4) |

where is the received pilot vector, is the pilot matrix, is the channel vector, is the activity vector, and refers to the element-wise multiplication. Therefore, this joint UAD and CE problem is formulated as follows

(5) |

where is the average propagation fading measured via statistical methods. According to assumptions (ii)-(iv), both the locations and transmission elevation angles of the stationary devices are fixed, which ensures the feasibility for the statistical measurement of and other static parameters , , and .

###### Remark 1

As in [7, 8], we employ Ku-band for the TSL in LEO satellite-enabled IoT. In the Ku-band, the propagation fading is related to the weather condition, the satellite receiver, and atmospheric variation [36]. We further assume that in each round of GF-RA, the LEO satellite only covers a relatively small area, where all devices experience similar weather condition due to geographical proximity. The assumption of geographical proximity can be justified by the fact that the coverage area of each LEO satellite will be greatly narrowed down in a forthcoming dense LEO constellation with thousands of satellites. In addition, each satellite can cover different small areas in a time-division manner, since IoT devices can be delay-tolerant in some typical applications [4, 13, 14]. In this way, in each round of GF-RA, all devices experience the same propagation fading due to geographical proximity. On the other hand, models the phase shift caused by orbital perturbation, non-ideal elimination of Doppler shift [40], or atmospheric variation [36]. Similarly, we can assume that all devices experience the same due to geographical proximity.

###### Remark 2

In order to reduce the access delay for this GF-RA system, we adopt non-orthogonal Gaussian pilot sequences. Specifically, we assume that each pilot symbol is independently drawn from a complex Gaussian distribution, i.e., . In contrast to orthogonal pilots with length , the length of non-orthogonal Gaussian pilots can be chosen flexibly. Furthermore, we employ pilots with length to reduce the access delay and improve the transmission efficiency of short data packets.

## Iii BR-MP-EM Algorithm for Joint UAD and CE

In order to address the joint UAD and CE problem for this terrestrial-satellite GF-RA scheme, we derive the BR-MP-EM algorithm in this section. The message passing process of this BR-MP-EM algorithm is illustrated on the factor graph in Fig. 2, where the -th sum node (SN) represents the summation term in (3), the -th variable node (VN) represents the product of the Bernoulli node and the Rician node for device , and the hyper-node represents the channel impairment term .

As shown in (2), the terrestrial-satellite channel of each device is a multiplication product of two terms, i.e. the Rician channel term and the channel impairment term . Accordingly, the message passing process on the factor graph is divided into two stages, i.e. the inner iterations (on the left sub-graph of Fig. 2) and the outer iterations (on the right sub-graph of Fig. 2). In the inner iterations, the propagation fading and phase shift are considered as known hyper-parameters, provided by the outer iterations. Then, we derive the Bernoulli-Rician message passing for the joint UAD and CE problem. In the outer iterations, the UAD and CE results from the left sub-graph are employed to update the estimates of hyper-parameters and with the EM method. Then another round of inner iterations is triggered by the updated hyper-parameters. The details of this BR-MP-EM algorithm are explained as follows.

### Iii-a Bernoulli-Rician Message Passing

In this subsection, we assume that the hyper-parameters are fixed as and after the -th outer iteration. As noted in [41] and in (2), the Rician channel of each device contains a LoS component and a complex Gaussian scattering component. Furthermore, as discussed in Section II-A, the LoS component is assumed deterministic for our proposed terrestrial-satellite GF-RA system. Therefore, the Rician channel follows the complex Gaussian distribution

(6) |

In this way, the Rician message of the channel estimation can be characterized by the estimation mean-value and estimation variance . In other words, is the estimate of in (6), while represents the deviation of this estimation. On the other hand, the Bernoulli message [34] for the activity of device is represented by a probability . That is, is the probability for to take the value 1. Then, we can derive the Bernoulli-Rician message passing for the joint UAD and CE problem on the left sub-graph of Fig. 2. The message passing diagram is illustrated in Fig. 3.

#### Iii-A1 Message Update at Sum Nodes

Denote as the Bernoulli message for the activity of device , which is passed from VN to SN in the -th inner iteration. Denote and as the Rician messages passed from VN to SN . That is, and represent the estimate and estimation deviation of the channel , respectively. Then we rewrite (3) as follows,

(7) |

where is approximated as an equivalent noise . Based on the incoming Bernoulli messages and Rician messages, we can derive and as follows,

(8) |

where represents the noise variance of in (7). represents the probability for the Bernoulli variable to take the value , i.e., . Then the messages passed from SN to VN are derived as follows.

Rician Message Update at SN: The Rician messages of , i.e. the estimate and estimation deviation passed from SN to VN in the -th inner iteration are derived as follows

(9) |

where and represent the expectation and variance of conditioned on , respectively.

Bernoulli Message Update at SN: The Bernoulli message of , i.e. the non-zero probability passed from SN to VN in the -th inner iteration is derived as follows

(10) |

where and represent the mean-value and variance of when , and

represents the probability density function (pdf) of a complex Gaussian distribution

, i.e.,(11) |

In order to avoid the computation overflow caused by a large number of multiplications of probabilities, we employ the log-likelihood ratio (LLR) to represent the Bernoulli message. The relationship between the LLR message and the non-zero probability is defined as

(12) |

where equation () of (12) is derived by substitution of (10).

#### Iii-A2 Message Update at Variable Nodes

The messages passed from VNs to SNs follow the message combination rule [38]. The Rician messages and Bernoulli messages are updated as follows.

Rician Message Update at VN: For the Rician messages, the estimate and estimation deviation of passed from VN to SN are derived as follows,

(13) |

where and represent the prior mean-value and variance of after the -th outer iteration, is the set of all the SNs, , and .

Bernoulli Message Update at VN: The Bernoulli message passed from VN to SN is derived as follows,

(14) |

where , and is the prior activation probability of each device. Again, for computational convenience, we derive the LLR message for the Bernoulli variable ,

(15) |

where is the prior LLR for the activity of each device. Note that, in the calculation of and in (8), the non-zero probability in (14) will be used for the -th SN update, and we have

(16) |

#### Iii-A3 CE Output and UAD Decision at Variable Nodes

After a fixed number (denoted by ) of inner iterations, each VN outputs its final CE and UAD decision. Similar to the message update at VNs, the CE output and UAD decision are obtained by the message combination rule over all the incoming messages.

CE Output: The final channel estimate and estimation deviation for device are obtained by combining all the incoming Rician messages about , i.e.,

(17) |

UAD Decision: The LLR for UAD decision is expressed as

(18) |

If , the activity is detected as , and then in (17) is considered as estimated channel for device . Otherwise, . Note that compared with (15), an additional LLR term is included in (18) to exploit the CE output to improve the UAD accuracy. The motivation to include and its derivation are detailed as follows.

Derivation of : Note that if the final estimate in (17) is close to , it is more likely that this device is inactive. Therefore, we can exploit the CE result in (17) to improve the UAD accuracy. The estimated channel can be expressed as , where represents the estimation error caused by channel noise and multi-user interference. Since in (17) represents the CE deviation, we approximate

as a complex Gaussian random variable with variance

, i.e. . Note that if , the prior information is that , and we further have . On the other hand, is equivalently equal to if , and then we have . Therefore, is derived as follows(19) |

where is defined in (11).

### Iii-B Expectation-Maximization Update of Hyper-Parameters

If device is detected as active after a certain number of inner iterations, and are passed to the right sub-graph of Fig. 2 to update the hyper-parameters and with the EM method. Note that the Rician messages and from the left sub-graph of Fig. 2 characterize the posterior complex Gaussian pdf of . For notational simplicity, we rewrite this posterior Gaussian pdf as ,

(20) |

It is shown that is obtained with given hyper-parameters and after the -th outer iteration. According to (6), we define the prior Gaussian pdf conditioned on hyper-parameters and ,

(21) |

Then, the EM update in the ()-th outer iteration aims to find a pair of hyper-parameters and that maximizes the expectation over the distribution , i.e.,

(22) |

where equation () of (22) is obtained by the fact that is a monotonic function, and equation () of (22) is obtained by the fact that the expectation is taken over with posterior pdf .

Considering that all the active devices share the same hyper-parameters, we further extend the EM update equation (22) from the single-user case to the multi-user case, i.e.,

(23) |

where is the set of devices that are detected as active and participate in the EM update. To solve the maximization problem in (23), we differentiate with respect to and respectively, and set both partial derivatives to . In this way, we have the following set of equations,

(24) |

By solving (24), we have the solution to (23). Then, the hyper-parameters in the ()-th outer iteration are updated as follows

(25) |

where represents the angle of a complex number , , . represents the averaging operation, i.e., where is the number of elements in the set . The detailed derivation of (25) is explained in Appendix A.

### Iii-C Summary of BR-MP-EM Algorithm

The proposed BR-MP-EM algorithm is summarized in Algorithm 1. In Algorithm 1, is the number of outer iterations, and is the number of inner iterations. That is, a total number of iterations are performed by the BR-MP-EM algorithm. Specifically, Line 3-4 represent the Rician message update and Bernoulli message update at SNs, respectively. Line 5-6 represent the Rician message update and Bernoulli message update at VNs, respectively. Line 8-9 represent the CE output and UAD decision of the inner iterations, respectively. Note that in Line 10, we calculate the relative variation to evaluate the convergence of the channel estimation. In line 11, we choose the set of devices for the EM update in the outer iteration. That is, the channel estimate of device will be included in the EM update if device is detected as active and the relative variation is smaller than a threshold . It is noted that the threshold is considered to avoid error propagation from inner iterations to outer iterations, and the value of () is chosen empirically. Experiments have shown that under different simulation configurations, setting could achieve better trade-off between convergence speed and robustness against error propagation. Line 12 represents the EM update in the outer iteration. In addition, Line 13-14 represent the message initialization for the next round of inner iterations after the -th EM update. Finally, the BR-MP-EM algorithm outputs the UAD result and CE result .

###### Remark 3

The proposed BR-MP-EM algorithm treats the random phase shift and propagation fading as hyper-parameters, and updates these hyper-parameters with the EM method. It is worth mentioning that the EM method does not rely on specific distributions of and . Instead, only the average measured propagation fading is needed for this algorithm.

###### Remark 4

In (2) we assume the channel impairments affect both the scattering component and the LoS component. However, this assumption does not sacrifice any loss of generality for the BR-MP-EM algorithm. That is, if the channel impairments only affect the LoS component, only some minor modifications are required on the BR-MP-EM algorithm in Algorithm 1. The details are provided in Appendix B.

###### Remark 5

It is assumed in Algorithm 1 that devices experience the same channel impairment . However, the BR-MP-EM algorithm can be readily generalized to the case when devices experience different channel impairments. Assume that devices are divided into groups by geographical proximity, and only the devices in the same group experience the same channel impairment. Since IoT devices are assumed stationary, the group membership can be acquired by the satellite. Then, in order to generalize the BR-MP-EM algorithm, we only need to independently perform the EM update (Line 12 of Algorithm 1) for each group of devices.

### Iii-D Computational Complexity

One remarkable property of the message passing algorithms is that the overall processing can be decomposed into many local computations, which are executed in parallel at different SNs and VNs on the factor graph. Therefore, the proposed BR-MP-EM algorithm exhibits low computational complexity, which is favorable for reducing the processing delay at the satellite receiver. The computational complexity of the proposed BR-MP-EM algorithm is analyzed as follows.

We evaluate the computational complexity by the number of real-number multiplication (or division) operations and the number of exponential (or logarithm) operations. In each inner iteration and outer iteration, the complexity for calculating related messages is listed in Table. I. It is shown that the computational complexity is mainly incurred by the Bernoulli-Rician message passing in the inner iterations, while the operations in the EM update are almost negligible. According to Table I, each inner iteration costs about multiplications (or divisions) and exponential (or logarithm) operations. Furthermore, as in [23], if we approximate the extrinsic messages ,