3D On and Off-Grid Dynamic Channel Tracking for Multiple UAVs and Satellite Communications

The space-air-ground integrated network (SAGIN) has drawn increasing attention for its benefits, such as wide coverage, high throughput for 5G and 6G communications. As one of the links, space-air communications between multiple unmanned aerial vehicles (UAVs) and Ka-band orbiting low earth orbit (LEO) satellites face a crucial challenge in tracking the 3D dynamic channel information. This paper exploits a statistical dynamic channel model called the multi-dimensional Markov model (MD-MM), which investigates the more realistic spatial and temporal correlation in the sparse UAVs-satellite channel. Specifically, the spatial and temporal probabilistic relationships of multi-user (MU) hidden support vector, single-user (SU) joint hidden support vector, and SU hidden value vector are investigated. The specific transition probabilities that connect the SU and MU hidden support vector for both azimuth and elevation directions are defined. Moreover, based on the proposed MD-MM, we derive a novel multi-dimensional dynamic turbo approximate message passing (MD-DTAMP) algorithm for tracking the 3D dynamic channel in multiple UAVs systems. Furthermore, we also develop a gradient update scheme to recursively find the azimuth and elevation offset for 3D off-grid estimation. Numerical results verify that the proposed algorithm shows superior 3D channel tracking performance with smaller pilot overhead and comparable complexity.


I. INTRODUCTION
With increasing communication traffic demands for emerging applications, such as the internet of things (IoT), cloud computing, autonomous cars and so on, future networks have been expected to provide more robust communication services. Hence, the space-air-ground integrated network (SAGIN) has been widely studied for its advantages such as large coverage, high throughput and strong resilience for 5G and beyond systems [1]. SAGIN has mainly three segments: space segment consists of satellite network, air segment consists of aerial network, and ground segment consists of terrestrial network [2]. Acting as the middle segment, aerial network with unmanned aerial vehicles (UAVs) is cost-effective. It can provide J. Yu is with the Department of Electronic Engineering and Computer Science, Queen Mary University of London, London E1 4NS, UK (email: jiadong.yu@qmul.ac.uk).
Y. Gao (corresponding author) is with the Department of Electrical and Electronic Engineering, University of Surrey, Surrey GU2 7XH, UK (e-mail: yue.gao@ieee.org).
X. Shen is with the Department of Electrical and Computer Engineering, University of Waterloo, Waterloo, ON N2L 3G1, Canada (e-mail: sshen@uwaterloo.ca). a shorter line of sight (LoS) range for transmission compared with satellite-ground links [3]. UAVs are deployed as base stations in both crowded areas and non-accessible emergency scenarios for UAV-aided air-ground communications, which can provide seamless coverage. For space-air links, UAVs play a role as a relay to transmit the information received from satellite (e.g. low earth orbit (LEO) and medium earth orbit (MEO) satellites [4]) to terrestrial network with mmWave Kaband [5], [6]. In the above two UAV-related communication links, due to the velocity of UAVs, it is undeniable that the estimation of the highly dynamic communication channels with lower overhead and good accuracy faces some challenges: the spatial and temporal variations caused by the mobility of UAVs in three dimensional (3D) space [7].

A. Related Work
Conventionally, channel state information (CSI) is estimated using techniques, such as least square (LS) and minimum mean square error (MMSE). However, for these methods, training overhead is proportional to the number of transmit antennas which is infeasible to mmWave communication with massive MIMO system [8]. To reduce the training overhead, compressive sensing (CS) based techniques, such as orthogonal matching pursuit (OMP) and compressive sampling matching pursuit (CoSaMP), have been widely studied to estimate the approximate sparse wireless channel in different domain [9].
Based on fundamental CS concepts and algorithms, some studies explored the more detailed sparse structure of the channel through spatial correlation [10]- [13] and temporal correlation [14]- [17]. Because of the physical scattering of the environment, [10] has concluded that the angular domain shows a burst sparsity structure. Hence, they proposed the burst least absolute shrinkage and selection operator (burst-LASSO) to extract the burst structure properties in the massive MIMO channels with enhanced performance. [11] also exploited the clustered sparse structure of the channel support in the spatial domain through Markov prior. [12] proposed a hidden Markov model to capture the angular domain structured sparsity of the channel. Moreover, [13] exploited the double sparsity of the multipath components in both angular domain and delay domain. However, these studies only considered the angular domain structure ( [13] also considered the delay domain) based on static channel environment. Effectively utilizing the prior information can reduce the training overhead and maintain a good current channel estimation performance for dynamic channel tracking. [14] proposed a prior aided channel tracking scheme by taking the practical user motion model into account through a temporal variation law of the physical direction. [15] modeled the channel by a Markov process and the Kalman filter to track the dynamic model. A structured matching pursuit [16] and differential orthogonal matching pursuit (D-OMP) [17] are proposed both by considering the temporal correlation that path delays change slowly over time but path gains evolve faster.
Jointly considering both spatial correlation and temporal correlation can improve the channel tracking performance [18]- [21]. [18] proposed the channel tracking scheme based on the prior-aided iterative angle and Doppler shift estimation algorithm and with the help of the previously estimated channel parameters. [19] decoupled the channel information into angles of arrival/departure (AoAs/AoDs) information and gain knowledge. A modified unscented Kalman filter tracks the AoAs/ AoDs, and the gain is estimated from the proposed beam training. [20] and [21] formulated the more realistic dynamic models of the AoDs support over time through hidden Markov model in flat and 3D space, respectively.

B. Motivation and Contribution
The aforementioned channel estimation and tracking work in [10]- [12], [14]- [20] can be directly adapted to the ground links and air-ground links with mobile scenario. However, for spaceair links, which consist of multiple UAVs and one satellite, there are some non-neglected features: the known orbiting satellite, the 3D UAVs moving trajectories, and the dominant LoS ray between each UAV and satellite [2].
Although [10]- [12] considered the angular domain sparsity of the channel, it's non-trivial to include the temporal information for channel tracking, especially in highly mobile space-air links with orbiting satellite and UAVs navigation. Comparatively, [14]- [17] all considered the prior channel information for time-varying channel tracking. However, they ignored the spatial angular structured sparsity. [18] studied the prior information of both angle and Doppler shifts for channel tracking, the probabilistic relationship of the angle movement is not explored. Moreover, although [19] and [20] both jointly considered spatial and temporal correlation of the dynamic channel, the premise of the system were both communicated through a flat surface. Furthermore, although our previous work [21] considered a more realistic 3D space-air communications with only one UAV, however, the multiple-UAVs system can provide wider coverage in practice. Compared with one UAV in the system, for multiple UAVs, there will have more than one structured sparse points on the quantized grids to represent the multi-user (MU) spatial correlation with sharing azimuth or elevation supports. Hence, a more complicated scenario with multiple UAVs should be further explored. The above algorithms all assumed that the true angles lie on the quantized grids. However, in the real world, the true angle is continuously distributed [22]. To solve the off-grid problem, [22] and [23] all proposed algorithms based on a two-dimensional (2D) flat environment. Hence, none of the mentioned work considered the 3D off-grid case of the channel model. Motivated by the above literature review, we consider the multiple UAVs space-air links to recursively track the 3D offgrid dynamic channel with smaller training overhead, comparable complexity, and improved estimation performance. The main contributions are summarized as follows: • We develop a statistical 3D on-grid and off-grid channel model called multi-dimensional Markov model (MD-MM) for space-air communication links with 3D navigation multiple UAVs and orbit LEO satellite in the system. The MD-MM captures the structured sparsity of the MU azimuth and elevation spatial domains and the probabilistic relationship between MU and SU support. • Based on the proposed MD-MM channel model, we propose a novel multi-dimensional dynamic turbo approximate message passing (MD-DTAMP) algorithm which can recursively track the 3D on-grid dynamic channel over time. • We further develop a gradient azimuth and elevation 3D off-grid estimation scheme for a more realistic channel model, which iteratively update the optimum 3D offset through gradient calculation. • The proposed MD-DTAMP algorithm with and without an off-grid scheme all show superior estimation performance compared with benchmark algorithms. Moreover, the proposed algorithm can achieve solid performance with smaller pilot overhead and comparable complexity. This paper is organized as follows. The downlink LEO satellite to multiple UAVs system model is introduced in Section II. In Section III, the 3D on-grid and off-grid MD-MM channel model are elaborated. The problem of channel tracking is formulated in Section IV. The proposed algorithm MD-DTAMP and the 3D gradient off-grid estimation scheme are presented in Section V. Numerical results and conclusions are then provided in Section VI and Section VII, respectively.
Notation: Throughout the paper, we use the following common notation. The complex numbers are denoted by C. The transpose and the conjugate transpose are denoted by (·) T , (·) H respectively. I is the identity matrix. CN (mean, covariance) indicates a complex Gaussian random vector with defined mean and covariance. · is the l 2 norm.

II. SYSTEM MODEL
In this paper, we consider the narrow-band flat fading transmission model through space-air downlink communication from LEO satellite to multiple UAVs. Since Ka-band is weather dependent, to mitigate the high path loss, the LEO satellite is installed with uniform planar array (UPA) with M = N x × N y antenna elements to compensate the transmission loss [6], [24]. The k th ∈ {1, ..., K} UAV is equipped with a high gain omnidirectional single antenna 1 [25], [26]. The illustration of the 1 For simplicity of only exploring the temporal and spatial correlation of the azimuth and elevation AoD of the satellite, the high gain omni-directional single antenna is assumed to be carried by the k th UAV. However, the proposed system model can be further extended with UAVs also installed with UPA. The main difference is that the channel model should further take the receiver array response vectors with azimuth and elevation directions of AoAs into consideration. This article has been accepted for publication in IEEE Transactions on Wireless Communications. This is the author's version which has not been fully edited and content may change prior to final publication. Citation information: DOI 10.1109/TWC.2021.3121301   I  TABLE OF IMPORTANT NOTATIONS   Notations Explanations n x , N x Index and total number of x-axis antenna/sparsity size n y , N y Index and total number of y-axis antenna/sparsity size N S Total pilot number m, M Index and total number of the antenna and the sparse grid t, T Index and total number of the time slot θ k,t , φ k,t k th UAV azimuth and elevation directions of AoD at t θ k,t ,φ k,t k th UAV on-grid azimuth and elevation directions of AoD at t ∆θ k,t , ∆φ k,t k th UAV azimuth and elevation offset directions of AoD at t θ range The range of azimuth directions of AoD φ range The range of elevation directions of AoD a (θ k,t , φ k,t ) Transmit array response vector at t h k,t k th UAV frequency domain channel model at t S k k th UAV transmit pilot symbols y RX k,t , y t k th UAV received signal at t y (t) k k th UAV total received signal for total period t n RX k,t , n k,t k th UAV additive white Gaussian noise at t g k,t k th UAV complex channel gain at t g k,t,m k th UAV channel gain element for the m th AoD at t g k,t k th UAV sparse matrix complex channel gain at t g (t) k , g (T ) k k th UAV total sparse matrix complex channel gain for period t or T b k,t,m k th UAV hidden support element for the m th AoD at t b k,t k th UAV hidden support vector at t b k th UAV hidden support vector for total period t or T b t MU Hidden support vector at t a k,t,m k th UAV azimuth support element for the m th AoD at t e k,t,m k th UAV elevation support element for the m th AoD at t a k,t , e k,t k th UAV azimuth and elevation support vector at t a (t) k , a (T ) k k th UAV azimuth support vector for total period t or T e (t) k , e (T ) k k th UAV elevation support vector for total period t or T a t,m ,ẽ t,m MU azimuth and elevation support element for the m th AoD at t a t ,ẽ t MU azimuth and elevation support vector at t a (t) ,ã (T ) MU azimuth support vector for total period t or T e (t) ,ẽ (T ) MU elevation support vector for total period t or T ϑ k,t,m k th UAV hidden value element for the m th AoD at t ϑ k,t k th UAV hidden value vector at t ϑ (t) k , ϑ (T ) k k th UAV hidden value vector for total period t or T system model can be found in Fig. 1. The important notations of the system and channel model are shown in Table I. In training procedure, as symbols S k ∈ C M ×N S are transmitted from satellite to UAV, the received signal at k th UAV at time slot t can be expressed as with h k,t ∈ C 1×M as the frequency domain channel vector, and n RX ∼ CN (0, σ 2 n I) is the additive white Gaussian noise.

A. Channel Model
The LoS ray is the dominant path in space-air communication links [27]. Hence, the h k,t denotes the frequency domain channel vector with one dominant LoS path, which is written as, whereg k,t is the complex channel gain, a (θ k,t , φ k,t ) ∈ C (Nx×Ny)×1 represents the transmit array response vectors, θ k,t and φ k,t as azimuth and elevation directions of AoD for k th UAV at time slot t. As for UPA installed LEO satellite, the steering vector can be written as where λ is the signal wavelength, 0 n x < N x and 0 n y < N y are the indices of an antenna element with N x and N y as antenna array sizes.

B. On-grid and Off-grid Virtual Channel Model
Conventionally, to explore the sparse property of the channel, AoDs are assumed to lie on the quantized grids for simplicity [28]. With discretized AoDs, the channel model is called on-grid channel. However, in practice, the true angle is continuously distributed. The channel with this kind of AoDs is named as off-grid channel [23]. Hence in this subsection, we present both on-gird and off-grid cases.
For on-grid case, the channel model can be further written into a sparse representation as where g k,t ∈ C 1×M is a sparse vector with one non-zero elements (which means only one dominant path for the k th UAV), M is the total number of spatial grids, the quantized AoD pairs are considered, i.e.,θ k,t = [0, θ range /N x , ..., θ range (N x −1)/N x ] andφ k,t = [0, φ range /N y , ..., φ range (N y − 1)/N y ], which means the virtual k th UAV's AoD pair can consist of a M uniform grid, and A θ k,t ,φ k,t ∈ C M ×M is the dictionary matrix. For off-grid case, we define ∆θ k,t,nx and ∆φ k,t,ny to represent the offset angles for both azimuth and elevation directions, respectively. Hence, the true AoD azimuth and elevation angle should be the combination of the nearest on-grid angle and the relative offset This article has been accepted for publication in IEEE Transactions on Wireless Communications. This is the author's version which has not been fully edited and content may change prior to final publication. This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TWC.2021.3121301, IEEE Transactions on Wireless Communications 4 angle. As an example, we demonstrate the k th UAV's AoD azimuth direction. The illustration can be found in Fig. 2. We assume that the true azimuth angle θ k,t has the nearest on-grid angleθ k,t,nx and the offset angle is ∆θ k,t,nx = θ k,t −θ k,t,nx . Similarly, the AoD elevation direction offset angle can be written as ∆φ k,t,nx = φ k,t −φ k,t,ny . Hence, the off-grid channel can be written as whereÃ θ k,t + ∆θ k,t ,φ k,t + ∆φ k,t = a θ k,t,1 + ∆θ k,t,1 ,φ k,t,1 + ∆φ k,t,1 , ..., a θ k,t,nx + ∆θ k,t,nx ,φ k,t,ny + ∆φ k,t,ny , ..., a θ k,t,Nx + ∆θ k,t,Nx ,φ k,t,Ny + ∆φ k,t,Ny is the off-grid array response vector dictionary matrix.
To be noticed, the AoD azimuth offset vector The otherwise case represents the non-active virtual on-grid offsets which can be observed in Fig. 2 with the notation from ∆θ k,t,1 to ∆θ k,t,N without the active path offset ∆θ k,t,nx . Similarly, the AoD elevation offset vector ∆φ k,t = ∆φ k,t,1 ..., ∆φ k,t,Ny T with ∆φ k,t,ny = φ k,t −φ k,t,ny , Active path, 0, Otherwise.
When ∆θ k,t = 0 and ∆φ k,t = 0, the off-grid channel can be treated as the on-grid channel in (4).

III. MULTI-DIMENSIONAL MARKOV MULTIUSER
CHANNEL MODEL In our model, communication links connected between LEO satellite and multiple UAVs are considered. Because of the LEO satellite height, a certain portion of the earth's surface is covered at any given time slot. This means a constellation of LEO satellites is necessary for full coverage of the earth [29]. We consider a condition that one LEO satellite communicates with the same group of multiple UAVs simultaneously before any communication link is handed over to the adjacent LEO satellite. This means the LEO satellite would orbit in a certain distance range, limiting the azimuth and elevation angle in fixed ranges. With assumptions above, we transfer the complex continuous relative displacement relationships between LEO satellite and UAVs into time-varying angular domain variance with transition probability similar to [19]- [21] (this is further elaborated in Section III-C).
In the MU scenario with multiple UAVs, we specially explore the relationship between the MU and SU hidden support vector. The MU azimuthã (T ) and elevationẽ (T ) support vector are both modeled as two dimensional Markov model (2D-MM) similarly in [21].ã (T ) andẽ (T ) represent the MU azimuth and elevation support vector over total time period T . During time slot t,ã t = [ã t,1 , ...,ã t,Nx ] ∈ {0, 1} Nx withã t,nx = a 1,t,m ∨ a 2,t,m ∨ ... ∨ a K,t,m , which means MU azimuth support is the union of the SU azimuth support of total K UAVs at time slot t. Similarly, for MU elevation

A. Probability Model for SU Channel
The dynamic angular domain channel for the k th UAV g (T ) k = {g k,1 , ..., g k,T } with g k,t = {g k,t,1 , ..., g k,t,M } can be modeled as a probabilistic signal model with two hidden random processes, SU joint hidden support vector b , 1} M represents the SU joint hidden support vector at time t, which indicates the k th UAV channel sparsity with b k,t,m ∈ {1, 0}. The SU hidden value vector ϑ k,t = [ϑ k,t,1 , ..., ϑ k,t,M ] ∈ C M represents the temporal correlation of the k th UAV channel coefficients. The dynamic channel element for the m th AoD direction in the time t of the k th UAV can be written as where b k,t,m denotes whether there is an active path or not (i.e., b k,t,m = 1 means the m th AoD path is activated) and ϑ k,t,m represents the path gain 2 .
To formulate the more realistic dynamic channel (i.e., the inactive AoD path at the current time t may be activated at 1536-1276 (c) 2021 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TWC.2021.3121301, IEEE Transactions on Wireless Communications 5 time t + 1), a probabilistic channel model of k th UAV with channel prior distribution can be formed as where the k th UAV channel vector conditional prior is with δ(·) is the Dirac delta function. The factor graph of the combined k th UAV channel can be found in Fig. 3, with π k,t,m denotes the conditional prior p (g k,t,m | b k,t,m , ϑ k,t,m ). This factor graph shows that the final k th UAV channel gain g k,t is the combination of the SU joint hidden support vector b k,t (on the left hand side of Fig. 3) and the SU hidden value vector ϑ k,t (on the right hand side of Fig. 3).

B. Gaussian Markov Model for SU Hidden Value Vector
The illustration of the factor graph for the k th UAV SU hidden value vector ϑ k,t and the Gaussian-Markov model can be found in Fig. 3 (on the right hand side). Due to the path gains change smoothly over time [17], the hidden value vector can be formulated as the Gauss-Markov processes [30] where β k ∈ [0, 1], ω k,t,m ∼ CN (0, ζ k ), µ k ∈ C is the mean of the process. When β k = 0, then ϑ k,t,m = ϑ k,t−1,m , which indicates that ϑ k,t remain unchanged over time. When β k = 1, then ϑ k,t,m = β k ω k,t,m + µ k , which means that the SU hidden value vector shows the i.i.d Gaussian distribution with mean µ k over time. If 0 < β k < 1, the conditional probability can be written as The distribution of the steady state of the process is where p (ϑ k,1,m ) and p (ϑ k,t,m | ϑ k,t−1,m ) are the transition probabilities that equal to the factor nodes f k,1,m and f k,t,m in Fig. 3.

C. 2D-MM for Azimuth and Elevation MU Hidden Support Vector
The details of the azimuth and elevation MU hidden support vector are introduced as follows: To explain the 2D-MM for MU hidden support vector of azimuth and elevation directions, we take three UAV users as an example in Fig. 4(a). The relationship of the MU azimuth supportã t and the MU elevation supportẽ t with three activated paths are shown. To be noticed, in this example, there are two users share the same MU azimuth support. Hence, the MU azimuth supportã t has only two activated points at time t. For MU elevation support e t has three activated points at time t.
Because of the navigation of the UAVs [31], together with the known orbital voyage of the LEO satellite, the relative displacement of the UAV and satellite is continous. Hence, the spatial sparsity pattern will act highly related with the prior adjacent pattern, i.e.,ã t+1,nx has higher probability be influenced by the adjacentã t,nx−1 in previous time slot [10]. [17] has also verified that the channel support change slowly over time, which indicatesã t+1,nx has higher probability to be influenced by itself in previous time slotã t,nx . This relationship which considers both spatial correlation and temporal correlation is named as dynamic 2D sparsity. The illustration of this 2D-MM can be found in the factor graph in Fig. 4(b). The model over total time period T can be expressed as: where

D. 3D-HMM for azimuth and elevation SU joint hidden support vector
To explore the relationship between MU and SU, we modeled the joint prior distribution of MU and SU azimuth or elevation support as HMM with conditional probability. Take k th UAV SU azimuth support as an example: where ς (A) k,ãt,n x ,a k,t,nx is the transition probability equal to the factor node k,t,nx in Fig.5(b). Furthermore,ã (T ) and a (T ) k are the MU azimuth support vector and k th UAV SU azimuth support vector over time period T , respectively. Although in practice, multiple UAVs may share partially MU azimuth or elevation support, i.e., in Fig.5, user 2 and user 3 both have the 4 th grid activated which means the same azimuth support is shared, our HMM model can still apply. To be noticed, each element in a (T ) k is generated by the MU azimuth support vector in the HMM with certain conditional probability, i.e., k,t,1,1 . In another word, when the MU azimuth support elementã t,nx = 0, there is no way that the k th UAV SU azimuth support element a k,t,nx = 1. However, when the MU azimuth support elementã t,nx = 1, there is a certain probability ς (A) k,t,1,1 that the n th x element a k,t,nx in the k th UAV SU azimuth support vector is 1. The k th UAV SU elevation support can be modeled similarly.
Hence, the 3D-HMM of the SU joint hidden support b (T ) k over total period time T which combines both azimuth and elevation direction supports can be found in Fig. 5(c). The joint conditional prior of the channel support probability of the k th UAV is given by Ny ny=1 p b k,t,nx,ny |a k,t,nx , e k,t,ny .
IV. 3D OFF-GRID DYNAMIC SPARSE CHANNEL TRACKING PROBLEM FORMULATION By taking the off-grid channel (5) into consideration, the received signal at k th UAV can be written as with y RX k,t ∈ C 1×N S . To convert the received signal into the standard CS model, we reshape g k,t from 1 × M to M × 1. Hence the received signal can be rewritten as with the sensing matrix Φ k,t ∆θ k,t , ∆φ k,t = S T k Ã H θ k,t + ∆θ k,t ,φ k,t + ∆φ k,t T ∈ C N S ×M , y k,t ∈ C N S ×1 . For the training phase, the transmitted pilot S k is designed based on the partial discrete Fourier transform (DFT) random permutation measurement matrix Φ k,t ∆θ k,t , ∆φ k,t [11], [32]. Details are shown in Appendix A. Given the received signals of total t time period y (t) k , the purpose is to track the dynamic channel vector g k,t at time t. The channel vector can be estimated aŝ g k,t,m = E g k,t,m |y (t) k , ∆θ k,t , ∆φ k,t , when the expectation is over the marginal posterior p g k,t,m |y  (ML) as follows, V. MULTI-DIMENSIONAL DYNAMIC TURBO APPROXIMATE MESSAGE PASSING ALGORITHM WITH 3D OFF-GRID ESTIMATION Approximate Message Passing (AMP) based algorithms have been widely studied for user detection [33], [34] and sparse channel estimation [30], [35], [36] because of its good reconstruction accuracy and low complexity. Traditionally, AMP is applied based on a static prior information. However, the prior distribution knowledge is time-varying in our structured dynamic channel model. As can be seen from the factor graph Fig. 6, the complex structure makes it more challenging to gain the exact marginal posterior in (20). We further employ a turbo iterative framework [37] to improve the reconstruction performance. The details of the proposed algorithm is elaborated in this section.
The factor graph of the total MD-MM distribution Fig. 6. This factor graph contains three sub-graphs: the MU azimuth and elevation vectorã t andẽ t (on the left-hand side of Fig. 6), the SU hidden support vector b k,t (in the middle of Fig. 6), and the k th UAV SU hidden value vector ϑ k,t (on the right-hand side of Fig.  6) all through time T . The MU azimuth and elevation vector a t andẽ t are both structured as 2D-MM. The k th UAV SU hidden support vector b k,t is structured as 3D-HMM which contains both SU azimuth and elevation vector a k,t and e k,t .
The details of all the necessary nodes in the factor graph are already explained in Section III. The expression of the factor nodes, distributions and functional forms are shown in Table  II.
Due to the complex structure of the total factor graph of the MD-MM, it is challenging to achieve the exact marginal posterior in (20). In this section, to explore the dynamic channel structure, we explain the structural details of the proposed algorithm based on the turbo framework with the proposed MD-MM MMSE denoiser. In [37], conventional turbo compressed sensing is achieved based on a two-step iteration framework : linear minimum mean-square error (LMMSE) estimator and i.i.d. prior-based MMSE estimator. Our proposed algorithm will be introduced in four parts: Module A LMMSE Estimator, Module B MD-MM MMSE denoiser, message passing through time, and 3D off-grid estimation.

A. Module A: LMMSE Estimator
In Module A, the angular channel vector g k,t at time t for k th UAV is estimated based on the the received signal y k,t with the prior distribution CN g k,t ; g pri A,k,t , v pri A,k,t I . The extrinsic mean and variance are defined as [37] g post A,k,t = g pri and This article has been accepted for publication in IEEE Transactions on Wireless Communications. This is the author's version which has not been fully edited and content may change prior to final publication.  p(g pri B,k,t,m |g k,t,m ) CN g k,t ; g pri B,k,t , v pri B,k,t I π k,t,m g k,t,m , b k,t,m , ϑ k,t,m p g k,t,m |b k,t,m , ϑ k,t,m δ g k,t,m − b k,t,m ϑ k,t,m f k,1,m ϑ k,1,m k,t,m (ãt,n x , a k,t,nx ) p(a k,t,nx |ãt,n x ) p(a k,t,nx = 1|ãt,n x = 0) = 0, p(a k,t,nx = 1|ãt,n x = 1) = ς (A) k,t,1,1 k,t,m (ẽt,n x , e k,t,nx ) p(e k,t,nx |ẽt,n x ) p(e k,t,nx = 1|ẽt,n x = 0) = 0, p(e k,t,nx = 1|ẽt,n respectively. The extrinsic distribution of g k,t is given by where the extrinsic mean and variance are respectively given by and

B. Module B: Message Passing MMSE Denoiser
Following LMMSE estimator in module A, the main task for module B is to exploit the k th UAV spatial sparsity details of the proposed MD-MM in the multi-UAVs system. Hence, the extrinsic mean and variance pass through the graph, which means g pri B,k,t = g ext A,k,t and v pri B,k,t = v ext A,k,t [37]. In module B, the purpose is to calculate the approximate posterior distributions of p(g k,t,m |g pri B,k ) instead of (20) using the sum-product message passing structure shown in the factor graph. For message passing algorithm, the basic assumption is that g pri B,k,t = g k,t + z k,t where z k,t ∼ CN (0, v pri B,k,t I) is independent of g k,t [37].
The formulation details are elaborated in Appendix B and the message passing procedure can be summarized as follows: • Firstly, in module A, based on the prior information passed from the last time slot t, the message passing over the k th UAV Gaussian hidden value vector ϑ k,t are given by (52) with µ dyna and (σ dyna ) 2 as the input to update the message υ ϑ k,t,m →π k,t,m . • Secondly, in module B, g pri B,k,t and v pri B,k,t from module A are passed over the path g k,t,m → π t,m → b k,t,m → u k,t,m → a k,t,nx using (53) to (56). • Thirdly, the forward-backward message passing is executed over the proposed MD-MM MU azimuth branch with three paths: a k,t,nx → k,t,nx →ã t,nx , the 2D-MM of the MU azimuth supportã t , andã t,nx → k,t,nx → a k,t,nx using (59) to (68). • Then the message pass through path a t,nx → u t,m → e t,ny (using (69) to (71)) to another MD-MM elevation branch with similar three paths: e k,t,ny → k,t,ny → e t,ny , the 2D-MM of the MU elevation supportẽ t , and e t,ny → k,t,ny → e k,t,ny using (71) to (82). The overall details of the MD-MM azimuth and elevation SU joint hidden support estimation can be found in Algorithm 1. • The message is finally passed back over the path e k,t,ny → u k,t,m → b k,t,m → π k,t,m → g k,t,m using (83) to (86).
This article has been accepted for publication in IEEE Transactions on Wireless Communications. This is the author's version which has not been fully edited and content may change prior to final publication. Based on the calculated updated messages υ π k,t,m →g k,t,m , the posterior distributions can be written as p g k,t,m |g pri B,k,t ∝ υ π k,t,m →g k,t,m (g k,t,m ) υ c k,t,m →g k,t,m (g k,t,m ) , (27) where υ c k,t,m →g k,t,m (g k,t,m ) = CN g k,t,m ; g pri B,k,t,m , v pri B,k,t . Then the posterior mean and variance of each element of g k,t can be updated as and where Var g k,t,m |g pri B,k,t denotes the conditional variance of g k,t,m given g pri B,k,t . Then similar to (25) and (26), the extrinsic update mean and covariance can be written as and

C. Message Passing Through Time
For dynamic scenario, prior information which is also treated as the temporal correlation can be considered for better reconstruction of the time-varying channel model. Hence, we further consider the message passing through time for channel tracking. There are mainly three branches passing through time in Fig. 6 the factor graph of the MD-MM: the path of the k th UAV Gaussian hidden value vector ϑ k,t , the path of MU azimuth supportã t and the path of MU elevation supportẽ t .
1) The k th UAV Gaussian Hidden Value Vector ϑ k,t Passing Through Time: The message passing from factor node π k,t,m to variable node ϑ k,t,m should be υ true π k,t,m →ϑ k,t,m (ϑ k,t,m ) = ρ out k,t,m CN ϑ k,t,m ; g pri B,k,t , v pri B,k,t To be noticed, if b k,t,m = 0, g k,t,m would be 0 in (9), which makes ϑ k,t,m unobservable. Hence, similarly in [30], a threshold which is slightly less than 1 is introduced. The modified message passing is given by  Initialization:γ f 12: for n x = 2, ..., N x do k,t,1,1 . #Message Passing Over k th SU Azimuth to Elevation Support: 20: Calculate individual message ρ out a k,t ,k,t,nx,ny using (70). 21: Calculate message υ a k,t,nx →u k,t,m (a k,t,nx ) and υ u k,t,m →e k,t,ny (e k,t,ny ) as in (69) and (71), respectively. 22: Calculate υ u k,t,m →e k,t,ny as in (73) and υ k,t,ny →ẽt,n y ∝ e k,t,ny as in (82). # HMM k th SU Elevation Support Estimation Based on 2D-MM MU Elevation: 23: Similar to the procedures from step 2 to step 18 by using elevation related parameters. Details can be found from (76) to (79). 24 , ρ out k,t,m > Th, (34) where is set as a small positive value close to 0. Then the message passing from factor node f k,t+1,m to variable node This article has been accepted for publication in IEEE Transactions on Wireless Communications. This is the author's version which has not been fully edited and content may change prior to final publication. Citation information: DOI 10.1109/TWC.2021.3121301 ϑ k,t+1,m is formed as where and 2) The MU Azimuth Supportã t Passing Through Time: The azimuth support prior information passing to the next time υã t,nx →dt+1,n x (ã t,nx ) can be given as with 3) The MU Elevation Supportẽ t Passing Through Time: Similarly, the elevation support prior information passing to the next time can be written as υẽ t,ny →rt+1,n y (ẽ t,ny ) = with . (41)

D. 3D Off-grid Estimation
As there is no closed-form expression due to the complex MD-MM structure of g (t) k , it's hard to directly maximize the log-likelihood function ln p y (t) k , ∆θ k,t , ∆φ k,t in (21). We extend the proposed in-exact majorization-minimization (MM) algorithm in [38] to help us solve the off-grid problem. As MM is the generalization of the expectation maximization method, the convergence of expectation maximization algorithm can   can be updated similarly with known (i + 1) th iteration azimuth offset ∆θ prove that the in-exact MM can converge to a stationary solution [39]. Firstly, we iteratively construct the lower bound for ln p y (t) k , ∆θ k,t , ∆φ k,t which can be named as surrogate function, i.e., for any fixed azimuth and elevation offset pairs ∆θ k,t , ∆φ k,t , the surrogate function can be defined as [38] L ∆θ k,t , ∆φ k,t |∆θ k,t , ∆φ k,t = p g k,t |y which satisfies the following properties: L ∆θ k,t , ∆φ k,t |∆θ k,t , ∆φ k,t ln p y (t) k , ∆θ k,t , ∆φ k,t , ∀∆θ k,t , ∆φ k,t , L ∆θ k,t , ∆φ k,t |∆θ k,t , ∆φ k,t = ln p y ∂L(∆θ k,t , ∆φ k,t |∆θ k,t , ∆φ k,t ) ∂∆θ k,t Then ∆θ k,t and ∆φ k,t can be updated at each iteration as This article has been accepted for publication in IEEE Transactions on Wireless Communications. This is the author's version which has not been fully edited and content may change prior to final publication. Citation information: DOI 10.1109/TWC.2021.3121301 with (·) (i) represents for the i (th) iteration. To be noticed, for the (i + 1) th iteration elevation update, we use the updated (i + 1) th azimuth offset as fixed known knowledge in (48). However, the problems in (47) and (48) are non-convex. Hence the gradient update can be applied as follows: and where ∆

VI. NUMERICAL RESULTS
In this section, the complexity and the estimation performance of the proposed algorithm MD-DTAMP with and without 3D off-grid estimation are evaluated. We assume that the proposed algorithm is performed in the distributed UAVs. For the space-air link of each single UAV, the sparsity ratio for both azimuth and elevation directions is Ny which indicates only one dominant LoS ray for each UAV-satellite link. The parameters are set as follows 3 : We assume that there are total 3 UAVs in the system. The MU azimuth and elevation spatial correlation parameters {q  1 3 , 0.5, 1, 0.9375, 0.9375}, Th = 1 − 10 −2 , = 10 −7 , γ b k,t = 1 , the azimuth range θ range = π and the elevation range φ range = π/4, 3D Off-grid azimuth and elevation update step sizes are ∆ To be noticed, in our simulation, our proposed algorithm MD-DTAMP with 3D off-grid estimation is presented in Algorithm 3. MD-DTAMP on grid which is Algorithm 3 without step 15. The benchmark algorithms which can be employed in our MU 3D channel model are sparse Bayesian learning 3 The parameters can be trained through expectation maximization algorithm which can refer to [30] for details. The MU azimuth and elevation related parameters for mutiple UAVs system in this paper are referred 2D settings in [20].

3:
while not converge do #Module A LMMSE estimator: 4: Update g post A,k,t and v post A,k,t in (22) and (23), respectively.

5:
Calculate g pri B,k,t = g ext A,k,t , v pri B,k,t = v ext A,k,t using (25) and (26). #Module B MD-MM MMSE denoiser: 6: Message passing over Gaussian hidden value vector ϑ k,t using (52). 7: Message Passing Over g k,t,m → π k,t,m → b k,t,m → u k,t,m → a k,t,nx using (53)-(58) 8: AoD SU azimuth support a k,t and SU elevation support e k,t estimation in Algorithm 1.

11:
Update g post B,k,t,m and v post B,k,t using (28) and (29). 12: Update g pri A,k,t = g ext B,k,t and v pri A,k,t = v ext B,k,t using (30) and (31). 13: Repeat Module A and Module B until convergence. 14: Hencep g k,t |y      Update dynamic messages to the next time slot (t + 1): υ a k,t,nx →d k,t+1,nx , υ e k,t,ny →r k,t+1,ny and υ f k,t+1,m →ϑ k,t+1,m using (38), (40) and (35). 18: end for (SBL) and D-OMP proposed in [38] and [17], respectively. All simulations are average of 500 times.   O N S M 2 per iteration 4 . Fig. 7 shows the computation time of various algorithms versus the number of antennas. It can be observed that the proposed on-gird MD-DTAMP shows comparable complexity compared with D-OMP. Furthermore, compared with SBL algorithm which requires matrix multiplication and inversion with high computation complexity, the proposed MD-DTAMP with off-grid estimation consumes less computation time for different number of antennas.

B. Simulation Results
We choose time-averaged normalized mean square error (TNMSE) as the performance metric as [30] whereĝ t is the estimated result of g t in time t. This performance metric is to evaluate the recovery performance of various algorithms. In Fig. 8, we further compare the TNMSE performance metric versus SNR among various algorithms with antenna array elements set as M = 64, 256. With SNR increases, all four algorithms show better recovery performance. To be noticed, our proposed MD-DTAMP on-grid and off-grid tracking algorithms always show superior performance compared to benchmark algorithms. Moreover, the off-grid can further improve the performance through gradient update of the 3D offset.
In Fig. 9, we further consider the relationship between the tracking performance and the training pilot overhead. It is easy to be noticed that our proposed algorithm with or without offgrid estimation can maintain a good performance with smaller pilot overhead. This proves that our proposed algorithms can achieve a solid performance by taking the dynamic channel spatial and temporal correlation into consideration.  1536-1276 (c) 2021 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. The proposed on-grid and off-grid algorithms are compared in Fig. 10. The off-grid estimation achieves sufficient performance gain over the on-grid algorithm. In Fig. 10 (a), the performance over continuous total time slot is compared. It is obvious that, with time increases, the performance show superior gain for both on-grid and off-grid algorithms. Moreover, it can be observed from Fig. 10 (a) and (b) that with M increases from 64 to 256, the tracking performance also improves. This is because with a larger antenna array and denser virtual gird, and the on-grid 3D tracking angle can be estimated with a narrower grid and closer to the true angle.

VII. CONCLUSION
This paper exploits the on-grid and off-grid 3D channel modelling and tracking for space-air communications with multiple UAVs and an LEO satellite in the system. The relative displacement between the orbit movement of the LEO satellite and the 3D trajectory of multiple UAVs was modelled based on the statistical dynamic channel model called MD-MM. The proposed MD-MM was composed of the insight probability relationship of the spatial and temporal correlation of the MU hidden support vector, SU joint hidden support vector, and SU hidden value vector. We further developed a novel MD-DTAMP algorithm to track the dynamic channel recursively. Furthermore, the gradient update for the 3D off-grid estimation scheme was proposed, which can iteratively find the azimuth and elevation offset. In addition, we compared the computation complexity between benchmark algorithms and the proposed algorithm with and without the off-grid scheme. Simulations verified that the proposed algorithm derived from the MD-MM dynamic channel model can better exploit the more realistic space-air 3D dynamic channel to achieve more accurate recovery performance over other benchmark algorithms with smaller pilot overhead.
passing is performed over the HMM a k,t for given input messages υ u k,t,nx ,ny →a k,t,nx .
We firstly introduce the details of the MU azimuth support estimation. The condition can be divided into two section, one is when t = 1, the other is when t > 1. To be notice, the message from factor node u t,m to variable node a t,nx is υ u k,t,m →a k,t,nx = ρ (A)in k,t,nx = ny ρ in a k,t ,k,t,nx,ny ny ρ in a k,t ,k,t,nx,ny + ny 1 − ρ in a k,t ,k,t,nx,ny , which is deriviated based on the the matched dimension M from the total gain, n x from the azimuth support and n y from the elevation support. The message from factor node k,t,nx to variable nodeã t,nx is υ k,t,nx →ãt,n x ∝ a k,t,nx nx υ u k,t,m →a k,t,nx (a k,t,nx )ς where When t = 1, the MU azimuth forward parameter is (64) The azimuth backward parameterγ (A)b t,nx = 1 1+γ (A) is given in (65). The final MU azimuth message which represents the nonzero probability can be written as υã t,nx → k,t,nx =ρ out at,k,nx δ(ã t,nx − 1) + (1 −ρ out at,k,nx )δ(ã t,nx ), whereρ out at,k,nx can be seen in (67). The message from factor node k,t,nx to variable node a k,t,nx is υ k,t,nx →a k,t,nx (a k,t,nx ) = ρ (A)out k,t,nx δ(a k,t,nx − 1) where ρ (A)out k,t,nx =ρ out at,k,nx ς (A) k,t,1,1 . 4) Message Passing Over the Path a k,t,nx → u k,t,m → e k,t,ny : The message from factor node a k,t,nx to variable node u k,t,m is υ a k,t,nx →u k,t,m (a k,t,nx ) = ρ out a k,t ,k,t,nx,ny δ(a k,t,nx − 1) + (1 − ρ out a k,t ,k,t,nx,ny )δ(a k,t,nx ), where ρ out a k,t ,k,t,nx,ny can be found in (70). The message from factor node u k,t,m to variable node e k,t,ny is υ u k,t,m →e k,t,ny (e k,t,ny ) = ρ in e k,t ,k,t,nx,ny δ e k,t,ny − 1 + 1 − ρ in e k,t ,k,t,nx,ny δ e k,t,ny , where ρ in e k,t ,k,t,nx,ny can be found in (72). 5) Message Passing Over e k,t,ny → k,t,ny →ẽ t,ny , the 2D-MM of the MU Elevation Supportẽ t,ny , andẽ t,ny → k,t,ny → e k,t,ny : The details can be derived similarly as the message passing process of the MU azimuth supportã t estimation. The message from factor node u k,t,m to variable node e k,t,ny is υ u k,t,m →e k,t,ny = ρ (E)in k,t,ny = nx ρ in e k,t ,k,t,nx,ny nx ρ in e k,t ,k,t,nx,ny + nx 1 − ρ in e k,t ,k,t,nx,ny , which is deriviated based on the the matched dimension M from the total gain, n x from the azimuth support and n y from the elevation support. The message from factor node k,t,ny to variable nodeẽ t,ny is υ k,t,ny →ẽt,n y ∝ e k,t,ny ny υ u k,t,m →e k,t,ny (e k,t,ny )ς (E) et,n y ,e k,t,ny where When t = 1, the MU elevation forward parameter is The final MU elevation message which represents the nonzero probability can be written as υẽ t,ny → k,t,ny =ρ out et,k,ny δ(ẽ t,ny − 1) + (1 −ρ out et,k,ny )δ(ẽ t,ny ), whereρ out et,k,ny can be seen in (81). The message from factor node k,t,ny to variable node e k,t,ny is υ k,t,ny →e k,t,ny (e k,t,ny ) = ρ (E)out k,t,ny δ(e k,t,ny − 1) + (1 − ρ (E)out k,t,ny )δ(e k,t,ny ), where ρ (E)out k,t,ny =ρ out et,k,ny ς (E) k,t,1,1 . 6) Message Passing Over the Path e k,t,ny → u k,t,m → b k,t,m → π k,t,m → g k,t,m : The message from factor node e k,t,ny to variable node u k,t,m is υ e k,t,ny →u k,t,m (e k,t,ny ) = ρ out e k,t ,k,t,nx,ny δ e k,t,ny − 1 + 1 − ρ out e k,t ,k,t,nx,ny δ e k,t,ny , where ρ out e k,t ,k,t,nx,ny can be found in (84). The message from factor node u k,t,m to variable node b k,t,m is υ u k,t,m →b k,t,m (b k,t,m ) = ρ out b k,t ,k,t,m δ(b k,t,m − 1) + 1 − ρ out b k,t ,k,t,m δ (b k,t,m ) , where ρ out b k,t ,k,t,m = ρ out a k,t ,k,t,nx,ny ρ out e k,t ,k,t,nx,ny γ b k,t . The message υ b k,t,m →π k,t,m from variable node b k,t,m to factor node π k,t,m is the same as υ u k,t,m →b k,t,m (b k,t,m ). The message from factor node π k,t,m back to variable node g k,t,m is υ π k,t,m →g k,t,m (g k,t,m ) = ρ out b k,t ,k,t,m CN g k,t,m ; µ out k,t,m , (σ out k,t,m ) 2 + 1 − ρ out b k,t ,k,t,m δ (g k,t,m ) .