Robust Channel Estimation Scheme for Multi-UAV MmWave MIMO Communication with Jittering

Abstract

In unmanned aerial vehicle (UAV)-assisted millimeter-wave (mmWave) communications, the communication performance is significantly degraded by UAV jitter. We formulate a UAV-assisted mmWave channel model with hybrid beamforming for the impacts of UAV jitter. Then, we derive the distribution of angle of arrivals/departures (AOAs/AODs) with random fluctuation of the UAV attitude angle. We develop an iterative reweight-based robust scheme as the super-resolution AOAs/AODs estimation method. Specifically, we introduce the partially adaptive momentum (Padam) estimation method to optimize the objective function of the jittering UAV mmWave massive multi-input multioutput (MIMO) system. Finally, compared with existing channel estimation schemes, the proposed UAV mmWave channel estimation method can achieve robust super-resolution performance in AOAs/AODs and path gains estimation with numerical results. Therefore, the proposed channel estimation scheme is very suitable for UAV mmWave massive MIMO communications with jittering.

1. Introduction

In commercial and civilian applications, unmanned aerial vehicles (UAVs) have received a lot of attention. UAVs will further spur the fast growth of different kinds of relational industries [1,2]. UAV communications also have better link qualities [3,4], wireless coverage [5], and high-throughput services for hotspot areas [6,7] than conventional cellular communications. The aerial radio access network (aerial-RAN) can be formed by UAVs. The lead UAV is responsible for communication with other UAVs (e.g., the followers) in the aerial-RAN [8,9]. UAVs work at high speeds in 3D space, along with low-latency and high-resilience constraints [10,11]. UAV communication can also support high data rates, which are required for future smart city applications. What is more, it will be significantly useful for performance improvements in future wireless sensor networks, such as reduced energy consumption, larger service coverage, and higher data transmissions [12].

In consideration of that, multi-UAV communication networks experience the severe challenge of increasing connections between leader UAVs and follower UAVs. It is the same as the connections between the base station and the dense-deployed user terminals in mobile communication networks. On one hand, there is huge demand for high-speed data transmission. On the other hand, there is a severe shortage of bandwidth resources [13]. As one of the key technologies to meet the increasing rate of demand for next-generation communication, millimeter-wave (mmWave) communication has received much attention in recent years [14]. By leveraging the abundant frequency spectrum resource in the mmWave frequency band, mmWave communication can support a higher data rate than that of conventional wireless communication systems. However, one bottleneck of the mmWave communication technology is that severe environmental conditions and hardware constraints must be considered, for instance, high path loss, reduced diffraction, implementation costs, and power consumption [15,16].

With the large number of antenna elements, the massive multi-input multioutput (MIMO) communication systems can significantly improve the spectral and energy efficiency in mmWave frequencies. Because of narrow and high-gain beams in small regions, the data rate and channel capacity are increased based on spatial multiplexing. Additionally, the massive MIMO mmWave communication system satisfies the low-latency requirement. What is more, in consideration of severe path loss, the beamforming realized by MIMO antenna arrays has been widely used. Beamforming techniques are the standard way to provide the necessary signal to noise ratio (SNR) gain, and they provide spatial multiplexing by using highly directional beams, especially in local coverage scenarios [17]. Only when the beams are properly aligned with the signal propagation paths is the beamforming gain maximized. The acquisition of accurate communication channel information at the base station (BS) is crucial for mmWave MIMO communications. As one kind of directional link, multi-UAV communications are realized by antenna arrays with a large number of elements, which are feasible for mmWave signaling [18]. Because of the small wavelength, it is possible to package a large number of antenna elements at both the leading UAVs side and the following UAVs side.

The perfect channel station information (CSI) of multi-UAV links is assumed in most existing works for UAV communications. However, in practical UAV communication systems, the communication platform is susceptible to strong wind gusts. Therefore, it introduces random body jittering in UAVs, such as, yaw jitter in the horizontal direction, or pitch jitter in the vertical direction [19,20]. The non-negligible channel estimation errors are introduced by body jitter, causing great energy consumption and data rate reduction. Taking the impact of UAV jitter into account, it is critically important to develop a channel estimation scheme for a multi-UAV mmWave massive MIMO system with hybrid analog–digital beamforming (HBF). Furthermore, the multi-UAV model employs novel channel estimation schemes based on HBF, improving data rates and quality of service (QOS).

In the HBF system, the number of antennas is larger than that of radio frequency (RF) chains; the digital baseband part cannot be directly connected to all the antennas. Thus, high-dimensional MIMO channel estimation is difficult to accurately obtain [21]. In the mmWave massive MIMO system with HBF, a compressive sensing scheme based on an angular channel sparsity grid is proposed in [22].

Exploiting the sparsity of mmWave channels, the sparse signal recovery problem is formulated for channel estimation. It is based on the parametric channel model with a quantized angle of arrivals/angle of departures (AoAs/AoDs), which is called an angle grid [23]. For channel estimation, the grid-based orthogonal matching pursuit (OMP) algorithm is proposed for estimating the channels of many common communication systems. Additionally, the grid-based OMP algorithm for estimating channels of hybrid MIMO systems is also developed [24]. The training overhead is reduced, but the AoAs/AoDs is discrete in angular domain, called “on-grid” AoAs/AoDs. In practice, AoAs/AoDs in reality is continuously distributed, called “off-grid” AoAs/AoDs. The AoAs/AoDs assumption based on discrete grid will lead to power leakage, which seriously reduces the accuracy of channel estimation and deteriorates channel recovery performance. To solve this resolution limitation caused by the on-grid angle estimation, the weighted iterative methods based on one-dimensional line spectrum estimation are proposed in [25]. The iterative reweighting (IR) method combining dictionary parameter learning and sparse signal recovery is proposed to estimate the off-grid AoAs/AoDs, extending to two-dimensional mmWave channel estimation [26].

The estimation of AoAs/AoDs is optimized by iteratively processing the surrogate function of the objective function. The weight parameters between sparsity and data error are optimized by the iterative updates [27]. In [28], the classical orthogonal matching pursuit (OMP) algorithm is used to recover the channel, transforming into the angular domain representation through standard spatial Fourier transform. The grid-based OMP algorithm for estimating channels of hybrid MIMO systems is modeled as the sparse signal recovery problem. In compressive sensing issues, the redundant dictionary is introduced, which consists of array response vectors with finely quantized angle grids. Additionally, the proposed angular resolution is much better than that of the virtual channel model; the coherence of the redundant dictionary is reduced by the nonuniformly distributed angle grids. What is more, the classical iterative algorithm, named approximate message passing (AMP), is also introduced to solve the channel estimation problem with low computational complexity, converting it to a sparse signal reconstruction problem [29].

In this paper, we address the above issues and summarize the main contributions of the work as follows.

• Super-resolution channel estimation model for multi-UAV jitter communication system: In practical UAV communication systems, multi-UAV transceivers flying in the air may encounter strong wind gusts, which leads to random body jittering with angular deviation. The estimation accuracy of the AoAs/AoDs between the leader UAVs and the follower UAVs, is influenced by jittering, which results in non-negligible AoAs/AoDs estimation errors. Because of the climatic conditions and electromagnetic interference, the information about the UAV locations may be imperfect [30]. Resulting from the UAVs’ location uncertainty, the additional path loss damages communication links between the leader UAVs and the follower UAVs. Thus, perfect CSI knowledge of the UAVs cannot be guarantee, and the system performance is degraded because of the imperfect CSI [7,31]. In view of this, multiple antennas performing beamforming can be employed to improve spectral efficiency in multi-UAV jitter communication systems. In this paper, the IR-based method is introduced to estimate the off-grid AoAs/AoDs, combining dictionary parameter learning and sparse signal recovery. Thus, the novel robust super-resolution channel estimation scheme for the multi-UAV jitter communication system moves the estimated AoAs/AoDs from the initial angular domain grid to the nongrid angular domain.
• Partially adaptive momentum (Padam) estimation method: The proposed IR-based algorithm is iterated from the grid domain, and every iteration searches new AoAs/AoDs estimations from the previous iteration, while the gradient descent algorithms are easy to converge to local optimal solutions. An algorithm that quickly converges to the optimal solution is in high demand. The Padam estimation method is introduced to improve the IR-based algorithm, optimizing the objective function of the jittering UAV communication system. Jointing the Adam/Amsgrad algorithm and the stochastic gradient descent (SGD) momentum, the Padam algorithm can well fill the generalization gap for adaptive gradient methods and control the adaptive level of the optimization procedure for the angle estimation scheme in this paper.
• Normalized mean square error (NMSE)-based robustness: In addition to commonly used performance indicators, we observe the jitter-robust performance of the proposed scheme based on NMSE-based robustness to prove which parameters influence beam alignment issues. Compared with the classical OMP and the AMP channel estimation schemes, the proposed IR-based channel estimation method with multi-UAV jitter can achieve more robust performance. Additionally, the jitter channel estimation accuracy is significantly improved.

In the remainder of this work, Section 2 describes the system model and problem description, including the multi-UAV mmWave massive MIMO system model based on HBF architecture and the formula of the channel estimation problem description for the jitter UAVs’ communication system. The proposed robust super-resolution multi-UAV MIMO channel estimation with jittering is investigated in Section 3, where the Padam estimation algorithm is given and the IR-based channel estimation scheme for multi-UAV mmWave massive MIMO system is optimized. Section 4 presents the numerical results for performance evaluation, showing that the improved scheme obviously increases the accuracy of the off-grid jitter AoAs/AoDs estimation. The robust performance of the proposed scheme is state-of-art, even in jitter channel estimation. The conclusions and our future work are in Section 5.

Notations: The boldface and lowercase letters denote vectors, while the boldface and capital letters denote matrices. The notation ${\parallel x\parallel }_F$ denotes the Frobenius norm of x, ${\parallel x\parallel }_0$ denotes the $l_0$-norm of x, and ${\parallel x\parallel }_2$ denotes the $l_2$-norm of x. The notation ${\left(·\right)}^T$ denotes the transpose of a matrix or a vector, ${\left(·\right)}^H$ denotes the Hermitian transpose of a matrix or a vector, and ${\left(·\right)}^{-1}$ denotes the matrix inversion of a matrix or a vector. $tr\left(·\right)$ denotes a trace matrix with elements in the vector, and $diag\left(·\right)$ denotes a diagonal matrix with elements in the vector. The notations ${\mathbb{C}}^{A\times B}$ and ${\mathbb{R}}^{A\times B}$, respectively, denote the space of $A\times B$ complex values matrices and real values matrices. The notation $\mathbb{E}\left[·\right]$ denotes the statistical expectation.

2. System Model and Problem Description

2.1. Multi-UAV MmWave Massive MIMO System Model

In the multi-UAV mmWave massive MIMO system, there is a lead UAV and U follower UAVs. For the ease of presentation, the MIMO communication link between the lead UAV with one follower UAV is taken as an example. The lead UAV is equipped with $N_T$ transmit antennas, and the follower UAVs are equipped with $N_R$ receiver antennas, respectively.

Due to the suitable beamforming, the beams between the lead UAV and the follower UAVs can be specified within the desired steering angle range. The received signal power at the follower UAVs is improved, and the interference caused by the undesired UAVs is reduced. By the hardware structure of the antenna array, the beamforming architectures in the transceivers are roughly divided into three categories, called hybrid beamforming, analog beamforming, and digital beamforming. In view of high power consumption, high cost, and high hardware complexity, the architecture of the fully digital beamforming is unsuited for the multi-UAV mmWave massive MIMO system, especially with a large number of antennas in the transceivers. While different from the digital beamforming architecture, only one RF chain is required in the analog beamforming architecture; phase shifters or switches in the analog domain are used to implement the analog beamforming architecture. Due to the analog beamforming architecture enabled by phase shifters, only the phase of the signal can be adjusted at each antenna element, so the degrees of freedom are limited. The hybrid beamforming architecture is proposed to combine the advantages of the digital and analog beamforming architectures. A reduced number of RF chains are used to reduce cost and energy consumption in the hybrid beamforming architecture. With the multistream transmission, the architecture can meet the overall performance requirements. Considering a reasonable tradeoff between the energy efficiency and the spectral efficiency of the multi-UAV mmWave massive MIMO system, the hybrid beamforming architecture is the best implementation solution, combing the advantages of the digital and analog beamforming architectures. In this paper, the hybrid analog–digital beamforming architecture is adopted, as shown in Figure 1, to address the hardware limitation challenge. For the UAV transceivers, the overall HBF architecture consists of a low-dimensional digital baseband precoder, high-dimensional analog RF chains implemented using simple analog components in the transmitter, a low-dimensional digital baseband combiner, and a high-dimensional analog RF chains in the receiver.

The lead UAV transmitter and the follower UAV receiver are, respectively, equipped with $N_T^{RF}$ and $N_R^{RF}$ RF chains, where $N_T^{RF}<N_T$ and $N_R^{RF}<N_R$. $\mathbf{x}$, as the $N_T^{RF}\times 1$ symbol vector, denotes the transmitted symbol with normalized average symbol energy. $E\left\{\mathbf{x}{\mathbf{x}}^H\right\}=I_{N_R}$. The receiver signal at the follower UAV can be given by

$$\mathbf{y}={\mathbf{Q}}^H\mathbf{HPx}+\mathbf{n},$$

(1)

where the receiver signal $\mathbf{y}\in {\mathbb{C}}^{N_R^{RF}\times 1}$, the hybrid combining matrix $\mathbf{Q}\in {\mathbb{C}}^{N_R^{RF}\times N_R}$, the channel matrix $\mathbf{H}\in {\mathbb{C}}^{N_R}\times N_T$, the HBF matrix $\mathbf{P}\in {\mathbb{C}}^{N_T}\times N_T^{RF}$, and the combined received noise $\mathbf{n}\in {\mathbb{C}}^{N_R^{RF}\times 1}$.

To reflect the sparse characteristics of the multi-UAV mmWave massive MIMO system channels, a typical geometric Saleh–Valenzuela channel model is widely adopted. The channel matrix can be given by

$$\mathbf{H}=\sum _{l=1}^L{\rho }_l{\mathbf{a}}_A\left({\varphi }_{\mathrm{A},l}^{},{\mathit{\theta }}_{\mathrm{A},l}^{}\right){\mathbf{a}}_D^H\left({\varphi }_{\mathrm{D},l}^{},{\mathit{\theta }}_{\mathrm{D},l}^{}\right).$$

(2)

The array response vectors ${\mathbf{a}}_A\left({\varphi }_{\mathrm{A},l}^{},{\mathit{\theta }}_{\mathrm{A},l}^{}\right)$ and ${\mathbf{a}}_D\left({\varphi }_{\mathrm{D},l}^{},{\mathit{\theta }}_{\mathrm{D},l}^{}\right)$ are determined by the complex path gain $z_l$. The random variables ${\varphi }_{\mathrm{A},l}^{}$ and ${\mathit{\theta }}_{\mathrm{A},l}^{}$, respectively, stand for the azimuth and elevation AOA for the l-th path. ${\varphi }_{\mathrm{D},l}^{}$ and ${\mathit{\theta }}_{\mathrm{D},l}^{}$ are the azimuth and elevation AOD for the l-th path. L is the number of channel propagation paths, and $L\ll min\left(N_R,N_T\right)$. ${\rho }_l$ is the complex path gain.

2.2. The Channel Estimation Problem Description for Jitter UAVs

To model the jitter UAVs under an observed line-of-sight (LoS) path, the time-varying angles ${\varphi }_{\mathrm{A},l}^{}$, ${\mathit{\theta }}_{\mathrm{A},l}^{}$, ${\varphi }_{\mathrm{D},l}^{}$, and ${\mathit{\theta }}_{\mathrm{D},l}^{}$ follow Gaussian random distribution. The azimuth AOD ${\varphi }_{\mathrm{D},l}^{}$, for instance, is given by

$${\varphi }_{\mathrm{D},l,t}^{}={\varphi }_{\mathrm{D},l,0}^{}+\sum _{i=1}^t{\mu }_i,$$

(3)

where ${\varphi }_{\mathrm{D},l,0}^{}\sim \mu \left(0,2\pi \right)$ is a randomly set initial angle. Additionally, ${\varphi }_{\mathrm{D},l,t}^{}$ follows a uniform distribution at time t. The jittering value ${\mu }_i\sim \mathcal{N}\left(0,{\sigma }_{\mu }^2\right)$ is a random angle increment and follows a normal distribution [32,33]. In practice, the jitter angle ${\mu }_i$ varies between ${10}^{-1}$ rad and ${10}^{-3}$ rad, and it depends on the climatic conditions and the UAV model, for instance, strong wind gusts and the specification of the wireless front-end. Assuming LOS between each leader UAV and follower UAV, the jitter channel model is shown in Figure 2. The red beam means the actual beam direction, which is impaired by a wind gust, and the blue beam means the desired beam direction. ${\varphi }_{\mathrm{D},l}^{}$ and ${\widehat{\varphi }}_{\mathrm{D},l}^{}$ denote the actual AOD and the estimation AOD, respectively, the same as the other three time-varying angles.

For the typical $N_1\times N_2^{}$ uniform planar arrays (UPAs), take ${\mathbf{a}}_D\left({\varphi }_{\mathrm{D},l}^{},{\mathit{\theta }}_{\mathrm{D},l}^{}\right)$ as an example.

$$\begin{array}{cc}\mathbf{a}\left({\varphi }_{D,l}^{},{\mathit{\theta }}_{D,l}^{}\right)=\hfill & {\left[1,e^{j2\pi dsin{\varphi }_{D,l}^{}sin{\mathit{\theta }}_{D,l}^{}/\lambda },...,e^{j2\pi \left(N_1-1\right)dsin{\varphi }_{D,l}^{}sin{\mathit{\theta }}_{D,l}^{}/\lambda }\right]}^T\\ & \otimes {\left[1,e^{j2\pi dcos{\mathit{\theta }}_{D,l}^{}/\lambda },...,e^{j2\pi \left(N_2-1\right)dcos{\mathit{\theta }}_{D,l}^{}/\lambda }\right]}^T\hfill \end{array}$$

(4)

where ⊗ stands for the Kronecker product, $\lambda$ means the wavelength, and d means the antenna spacing.

Defining $\varphi \stackrel{\Delta }{=}dsin\varphi sin\mathit{\theta }/\lambda$ and $\mathit{\theta }\stackrel{\Delta }{=}dcos\mathit{\theta }/\lambda$, $\mathbf{H}$ in Equation (2) is given by

$$\mathbf{H}={\mathbf{A}}_D\left({\mathit{\theta }}_D\right)diag\left(\rho \right){\mathbf{A}}_A^H\left({\mathit{\theta }}_A\right),$$

(5)

where ${\mathit{\theta }}_D={\left[{\varphi }_{D,1}^{},{\mathit{\theta }}_{D,1}^{},{\varphi }_{D,2}^{},{\mathit{\theta }}_{D,2}^{},...,{\varphi }_{D,L}^{},{\mathit{\theta }}_{D,L}^{}\right]}^T$, ${\mathit{\theta }}_A={\left[{\varphi }_{A,1}^{},{\mathit{\theta }}_{A,1}^{},{\varphi }_{A,2}^{},{\mathit{\theta }}_{A,2}^{},...,{\varphi }_{A,L}^{},{\mathit{\theta }}_{T,L}^{}\right]}^T$, $\rho ={\left[{\rho }_1,{\rho }_2,...,{\rho }_L\right]}^T$, ${\mathbf{A}}_D\left({\mathit{\theta }}_D\right)=\left[{\mathbf{a}}_D\left({\varphi }_{D,1}^{},{\mathit{\theta }}_{D,1}^{}\right){\mathbf{a}}_D\left({\varphi }_{D,2}^{},{\mathit{\theta }}_{D,2}^{}\right)\cdots {\mathbf{a}}_D\left({\varphi }_{D,L}^{},{\mathit{\theta }}_{D,L}^{}\right)\right]$, and ${\mathbf{A}}_A\left({\mathit{\theta }}_A\right)=\left[{\mathbf{a}}_A\left({\varphi }_{A,1}^{},{\mathit{\theta }}_{A,1}^{}\right){\mathbf{a}}_A\left({\varphi }_{A,2}^{},{\mathit{\theta }}_{A,2}^{}\right)\cdots {\mathbf{a}}_A\left({\varphi }_{A,L}^{},{\mathit{\theta }}_{A,L}^{}\right)\right]$.

For comparison, the steering vector of the uniform linear arrays (ULAs) can be given as

$$\begin{array}{c}\mathbf{a}\left(\varphi \right)={\left[1,e^{j2\pi dsin\varphi /\lambda },...,e^{j2\pi (N-1)dsin\varphi /\lambda }\right]}^T.\hfill \end{array}$$

(6)

Define $\mathbf{X}=\mathbf{Px}$, where the complex vector $\mathbf{x}$ has size $N_T\times 1$. It is assumed that $N_X$ pilot sequences are, respectively, sent by the antennas at the transmitting end, and $N_X<N_T^{}$. The pilot sequences are ${\mathbf{X}}_1,{\mathbf{X}}_2,\cdots ,{\mathbf{X}}_{N_X}$. ${\mathbf{y}}_p$ with $N_Y^{}$ dimensions is obtained by M time slots, and $N_Y^{}=M{N}_R^{RF}$. In the m-th time slot, the received pilot sequence ${\mathbf{y}}_{p,m}$ with $N_R^{RF}$ dimensions is obtained by the combining matrix ${\mathbf{W}}_m$. ${\mathbf{y}}_{p,m}$ and ${\mathbf{y}}_p$ can be, respectively, given by

$${\mathbf{y}}_{p,m}={\mathbf{W}}_m^H\mathbf{H}{\mathbf{X}}_p+{\mathbf{n}}_{p,m},$$

(7)

$${\mathbf{y}}_p={\mathbf{W}}_{}^H\mathbf{H}{\mathbf{X}}_p+{\mathbf{n}}_p,$$

(8)

where ${\mathbf{y}}_p={\left[{\mathbf{y}}_{p,1}^T,{\mathbf{y}}_{p,2}^T,\dots ,{\mathbf{y}}_{p,M}^T\right]}^T$ has size $N_Y\times 1$, noise ${\mathbf{n}}_p$ has size $N_Y\times 1$, and $\mathbf{W}=$$\left[{\mathbf{W}}_1,{\mathbf{W}}_2,\dots ,{\mathbf{W}}_{\mathrm{M}}\right]$ has size $N_R\times N_Y^{}$. Additionally, Equation (7) can be expressed in vector form as

$$\mathbf{Y}={\mathbf{W}}_{}^H\mathbf{HX}+\mathbf{N}.$$

(9)

The channel $\mathbf{H}$ estimation can be replaced by the estimation of ${\mathit{\theta }}_D$, ${\mathit{\theta }}_A$, path gains $\rho$, and the number of paths. Considering the sparsity of angle domain channel and unknown number of paths, the problem of jitter channel estimation can be deduced as

$$\underset{{\widehat{\mathit{\theta }}}_D,{\widehat{\mathit{\theta }}}_A,\widehat{\rho }}{min}{∥\widehat{\rho }∥}_0,s.t.{∥\mathit{Y}-{\mathit{W}}_{}^H\widehat{\mathit{H}}\mathit{X}∥}_{\mathrm{F}}\le \epsilon ,$$

(10)

where ${∥\widehat{\rho }∥}_0$ is equivalent to the path estimation $\widehat{L}$, $\epsilon$ is the error tolerance parameter, and $\widehat{\mathbf{H}}$ is the estimation of $\mathbf{H}$. Then, we convert the channel estimation into a constrained optimization problem as shown in Equation (10), and the constraint condition is the function of $\widehat{\rho }$, ${\widehat{\mathit{\theta }}}_D$, and ${\widehat{\mathit{\theta }}}_A$. Due to the loss of channel angular quantization, Equation (10) incurs serious estimation accuracy problems in the existing channel estimation schemes. We propose the channel estimation method based on the IR algorithm to improve channel estimation accuracy through the iterative optimization of AoAs/AoDs. ${\mathit{\theta }}_D$, ${\mathit{\theta }}_A$, and path gains $\rho$ in Equation (10) can be estimated by the improved method.

3. Proposed Robust Super-Resolution Channel Estimation with Jittering

We extend the IR-based method of one-dimensional spectral estimation to two-dimensional channel estimation by solving Equation (10). Additionally, the angle estimation is optimized by the Padam estimation method.

3.1. Padam Estimation Algorithm

In order to achieve better performance on these various stochastic optimization tasks, many adaptive gradient methods have been proposed. The Adagrad method with an adaptive learning rate for each individual dimension was one of the first of these proposed algorithms. It is used to motivate the research on adaptive gradient methods in the machine learning community. What is more, Adagrad can achieve a significant convergence increase for these sparse gradient situations, and even for nonsparse gradient settings [34].

Following the idea of adaptive learning rate, the RMSprop algorithm is improved by changing the arithmetic averages of the Adagrad algorithm to exponential averages. Although the RMSprop algorithm is an empirical method without theoretical verification, the exponential moving average variants of Adagrad are increased in the the RMSprop algorithm.

With simple manipulations, the formula of the exponential moving average variant of Adagrad is improved in the Adam algorithm, which combines both the idea of momentum acceleration and the RMSprop algorithm. For the sake of convergence analysis, the step size is also needed in the Adam algorithm, and the decaying step size and even constant step size all work well in the Adam algorithm. To simplify, the bias correction step size is ignored in this paper, since it does not have an impact on the argument. However, in some special circumstances, the Adam algorithm may exhibit nonconvergence issues. It cannot accumulate long-term memory of past gradients, increasing the effective learning rate in some cases. Then, the modified algorithm, called the Amsgrad algorithm, is adopted by adding a step to ensure the decay of the effective learning rate in the Adam algorithm.

In view of the Adagrad, Adam, and Amsgrad algorithms above, these theoretical guarantees on adaptive gradient schemes are just proved for convex functions. The Padam algorithm is improved for filling the generalization gap in the adaptive gradient methods, and the partial adaptive parameters $p\in \left(0,1/2\right]$ are introduced to control the adaptive level of the optimization procedure.

3.2. Optimization of the Channel Estimation Scheme with Jittering

In fact, the main difficulty in solving Equation (10) is that $l_0$-norm has low computational efficiency in finding the optimal solution. Then, we introduce a log-sum function to replace $l_0$-norm, efficiently finding the sparse solution of $\widehat{\rho }$. It can be given by

$$\underset{{\mathit{\theta }}_D,{\mathit{\theta }}_A,\rho }{min}\mathrm{F}\left(\rho \right)\stackrel{\Delta }{=}\sum _{l=0}^{L}log\left({\left|{\rho }_l\right|}^2+\delta \right),s.t.{∥\mathbf{Y}-{\mathbf{W}}_{}^H\widehat{\mathbf{H}}\mathbf{X}∥}_{\mathrm{F}}\le \epsilon ,$$

(11)

where $\delta >0$ and $\lambda >0$ are introduced to further deduce the problem (11) as an unconstrained optimization problem. It can be written as

$$\underset{{\mathit{\theta }}_D,{\mathit{\theta }}_A,\rho }{min}\mathrm{G}\left({\mathit{\theta }}_D,{\mathit{\theta }}_A,\rho \right)\stackrel{\Delta }{=}\sum _{l=0}^{L}log\left({\left|{\rho }_l\right|}^2+\delta \right)+\lambda {∥\mathbf{Y}-{\mathbf{W}}_{}^H\widehat{\mathbf{H}}\mathbf{X}∥}_{\mathrm{F}}^2.$$

(12)

Then, the iterative surrogate function is used to replace the log-sum function in (11). Its minimization is equivalent to that of $min\mathrm{G}\left({\mathit{\theta }}_D,{\mathit{\theta }}_A,\rho \right)$. Additionally, it can be given by

$$\underset{{\mathit{\theta }}_D,{\mathit{\theta }}_A,\rho }{min}{S}^{\left(i\right)}\left({\mathit{\theta }}_D,{\mathit{\theta }}_A,\rho \right)\stackrel{\Delta }{=}{\lambda }^{-1}{\rho }^H{\mathbf{Z}}^{\left(i\right)}\rho +{∥\mathbf{Y}-{\mathbf{W}}_{}^H\widehat{\mathbf{H}}\mathbf{X}∥}_{\mathrm{F}}^2.$$

(13)

${\mathit{Z}}^{\left(i\right)}$ is defined as

$${\mathit{Z}}^{\left(i\right)}\stackrel{\Delta }{=}diag\left(\frac{1}{{\left|{\widehat{\rho }}_1^{\left(i\right)}\right|}^2+\delta }\frac{1}{{\left|{\widehat{\rho }}_2^{\left(i\right)}\right|}^2+\delta }...\frac{1}{{\left|{\widehat{\rho }}_L^{\left(i\right)}\right|}^2+\delta }\right),$$

(14)

where ${\widehat{\rho }}_{}^{\left(i\right)}$ is the estimation of $\rho$ at the i -th iteration. Additionally, there are better estimations for ${\widehat{\mathit{\theta }}}_D^{\left(i+1\right)}$, ${\widehat{\mathit{\theta }}}_A^{\left(i+1\right)}$, and ${\widehat{\rho }}^{\left(i+1\right)}$ to make ${\mathbf{Z}}^{\left(i\right)}$ smaller. This means

$$S^{\left(i\right)}\left({\widehat{\mathit{\theta }}}_D^{\left(i+1\right)},{\widehat{\mathit{\theta }}}_A^{\left(i+1\right)},{\widehat{\rho }}^{\left(i+1\right)}\right)\le S^{\left(i\right)}\left({\widehat{\mathit{\theta }}}_D^{\left(i\right)},{\widehat{\mathit{\theta }}}_A^{\left(i\right)},{\widehat{\rho }}^{\left(i\right)}\right).$$

(15)

Then,

$$\begin{gathered} G\left({\widehat{\mathit{\theta }}}_D^{\left(i+1\right)},{\widehat{\mathit{\theta }}}_A^{\left(i+1\right)},{\widehat{\rho }}^{\left(i+1\right)}\right)-\lambda S^{\left(i\right)}\left({\widehat{\mathit{\theta }}}_D^{\left(i+1\right)},{\widehat{\mathit{\theta }}}_A^{\left(i+1\right)},{\widehat{\rho }}^{\left(i+1\right)}\right) \\ \le G\left({\widehat{\mathit{\theta }}}_D^{\left(i\right)},{\widehat{\mathit{\theta }}}_A^{\left(i\right)},{\widehat{\rho }}^{\left(i\right)}\right)-\lambda S^{\left(i\right)}\left({\widehat{\mathit{\theta }}}_D^{\left(i\right)},{\widehat{\mathit{\theta }}}_A^{\left(i\right)},{\widehat{\rho }}^{\left(i\right)}\right). \end{gathered}$$

(16)

In (16), $G-\lambda S^{\left(i\right)}$ reaches the maximum value when ${\widehat{\rho }}^{\left(i+1\right)}={\widehat{\rho }}^{\left(i\right)}$. Combining Equations (15) and (16), there is

$$G\left({\widehat{\mathit{\theta }}}_D^{\left(i+1\right)},{\widehat{\mathit{\theta }}}_A^{\left(i+1\right)},{\widehat{\rho }}^{\left(i+1\right)}\right)\le G\left({\widehat{\mathit{\theta }}}_D^{\left(i\right)},{\widehat{\mathit{\theta }}}_A^{\left(i\right)},{\widehat{\rho }}^{\left(i\right)}\right).$$

(17)

So, the minimization of G is equivalent to that of $S^{\left(i\right)}$. Optimizing $\rho$ in Equation (12), the optimal $\widehat{\rho }$ and $S^{\left(i\right)}$ can be given as

$$\begin{array}{c}{\rho }_{opt}^{\left(i\right)}\left({\mathit{\theta }}_D,{\mathit{\theta }}_A\right)\stackrel{\Delta }{=}arg\underset{\rho }{min}{S}^{\left(i\right)}\left({\mathit{\theta }}_D,{\mathit{\theta }}_A,\rho \right)\hfill \\ \begin{array}{ccc}\begin{array}{cc}& \end{array}& & \end{array}={\left({\lambda }^{-1}\mathbf{Z}+{\displaystyle \sum _{p=1}^{N_X}}{\mathit{K}}_p^H{\mathit{K}}_p\right)}^{-1}\left({\displaystyle \sum _{p=1}^{N_X}}{\mathit{K}}_p^H{\mathbf{y}}_p\right),\hfill \end{array}$$

(18)

$$\begin{array}{c}{S}_{opt}^{\left(i\right)}\left({\mathit{\theta }}_D,{\mathit{\theta }}_A\right)\stackrel{\Delta }{=}\underset{\rho }{min}{S}^{\left(i\right)}\left({\mathit{\theta }}_D,{\mathit{\theta }}_A,\rho \right)\hfill \\ \begin{array}{cc}\begin{array}{cc}\begin{array}{cc}& \end{array}& \end{array}& \begin{array}{c}=-{\left({\displaystyle \sum _{p=1}^{N_X}}{\mathit{K}}_p^H{\mathbf{y}}_p\right)}^H{\left({\lambda }^{-1}\mathbf{Z}+{\displaystyle \sum _{p=1}^{N_X}}{\mathit{K}}_p^H{\mathit{K}}_p\right)}^{-1}\hfill \\ \begin{array}{ccc}\begin{array}{cc}·\left({\displaystyle \sum _{p=1}^{N_X}}{\mathit{K}}_p^H{\mathbf{y}}_p\right)+{\displaystyle \sum _{p=1}^{N_X}}{\mathbf{y}}_p^H{\mathbf{y}}_p& \end{array}& & ,\end{array}\hfill \end{array}\end{array}\hfill \end{array}$$

(19)

where ${\mathbf{K}}_p={\mathbf{W}}^H{\mathbf{A}}_{\mathrm{R}}diag\left({\mathbf{A}}_{{}^{\mathrm{T}}}^H{\mathbf{x}}_p\right)$. Thus, problem (10) can be transformed into the optimization of normalized space angle ${\mathit{\theta }}_D,{\mathit{\theta }}_A$ in problems (18) and (19).

3.3. Padam-Based Robust Super-Resolution Channel Estimation Scheme

To solve problems (18) and (19) in the jittering UAV communication system, the proposed robust channel estimation scheme is described in Algorithm 1. In the IR-based method, $\lambda$ is updated in each iteration. A larger $\lambda$ leads to a smaller residue, which means a well-fitted estimation. However, a smaller $\lambda$ is poorly fitted, which means sparser estimation. Therefore, a larger $\lambda$ is chosen and updated by

$$\lambda =min\left(\alpha /e^{\left(i\right)},{\lambda }_{max}\right),$$

(20)

where $\alpha$ is constant scale factor, i is the number of iterations, and ${\lambda }_{max}$ is used to well control the problem, for example, $\lambda ={10}^{-8}$. As the squared residue in i-th iteration, $e^{\left(i\right)}$ can be given as

$$e^{\left(i\right)}={∥\mathit{Y}-{\mathit{W}}_{}^H{\mathit{A}}_D\left({\widehat{\mathit{\theta }}}_D^{\left(i\right)}\right)diag\left({\widehat{\rho }}^{\left(i\right)}\right){\mathit{A}}_A^H\left({\widehat{\mathit{\theta }}}_A^{\left(i\right)}\right)\mathit{X}∥}_{\mathrm{F}}^2.$$

(21)

The iteration is started at the angle domain grids in the proposed algorithm. Searching for new estimations ${\widehat{\mathit{\theta }}}_D^{\left(i+1\right)}$ and ${\widehat{\mathit{\theta }}}_A^{\left(i+1\right)}$ in the i-th iteration, let $S^{\left(i\right)}$ be smaller. The search task can be accomplished by Padam estimation method, and it is

$$\begin{array}{c}{\mathit{m}}^i=\alpha ·{\mathit{m}}^{\left(i-1\right)}+\left(1-\alpha \right)·{\nabla }_{{\mathit{\theta }}_D}{S}_{\mathrm{opt}}^{\left(i\right)}\left({\widehat{\mathit{\theta }}}_D^{\left(i\right)},{\widehat{\mathit{\theta }}}_A^{\left(i\right)}\right),\hfill \\ {\mathit{v}}^i=\beta ·{\mathit{v}}^{\left(i-1\right)}+\left(1-\beta \right)·{\left({\nabla }_{{\mathit{\theta }}_D}{S}_{\mathrm{opt}}^{\left(i\right)}\left({\widehat{\mathit{\theta }}}_D^{\left(i\right)},{\widehat{\mathit{\theta }}}_A^{\left(i\right)}\right)\right)}^2,\hfill \\ {\widehat{\mathit{v}}}^i=max\left({\widehat{\mathit{v}}}^{\left(i-1\right)},{\mathit{v}}^i\right),\hfill \end{array}$$

(22)

where ${\nabla }_{{\mathit{\theta }}_D}{S}_{\mathrm{opt}}^{\left(i\right)}\left({\widehat{\mathit{\theta }}}_D^{\left(i\right)},{\widehat{\mathit{\theta }}}_A^{\left(i\right)}\right)$ is the gradient, and $\alpha$ and $\beta$ are momentum parameters.

${\widehat{\mathit{\theta }}}_D^{(i+1)}$ and ${\widehat{\mathit{\theta }}}_A^{(i+1)}$ can be written as

$$\begin{gathered} {\widehat{\mathit{\theta }}}_D^{(i+1)}=\prod _{F,diag{\left({\widehat{v}}^i\right)}^p}\left({\widehat{\mathit{\theta }}}_D^{\left(i\right)}-\eta ·m^i/\phantom{m^i{\left({\widehat{v}}^i\right)}^p}{\left({\widehat{v}}^i\right)}^p\right) \\ {\widehat{\mathit{\theta }}}_A^{(i+1)}=\prod _{F,diag{\left({\widehat{v}}^i\right)}^p}\left({\widehat{\mathit{\theta }}}_A^{\left(i\right)}-\eta ·m^i/\phantom{m^i{\left({\widehat{v}}^i\right)}^p}{\left({\widehat{v}}^i\right)}^p\right) \end{gathered}$$

(23)

where $\eta$ is step length, F is the initial grid point, and $p\in \left(0,1/2\right]$ is the partial parameters.

Algorithm 1:The Proposed Robust Super-Resolution UAV Jitter Channel Estimation Scheme

Require:  Transmit signals $\mathbf{X}$, receive signals $\mathbf{Y}$, beamforming matrix $\mathbf{W}$, termination threshold ${\epsilon }_{th}$, initial ${\widehat{\mathit{\theta }}}_D^{\left(0\right)}$, ${\widehat{\mathit{\theta }}}_A^{\left(0\right)}$, and the threshold ${\rho }_{th}$ of path gains. Ensure:  The estimation of AoAs ${\widehat{\mathit{\theta }}}_A^{}$, AoDs ${\widehat{\mathit{\theta }}}_D^{}$, and path gains $\widehat{\rho }$ for all paths. 1: Set ${\widehat{\rho }}^{\left(0\right)}={\rho }_{opt}\left({\widehat{\mathit{\theta }}}_D^{\left(0\right)},{\widehat{\mathit{\theta }}}_A^{\left(0\right)}\right)$ by Equation (18). 2: repeat 3:     $\lambda$ is updated by Equation (20). 4:     According to Equation (19), set $S_{opt}^{\left(i\right)}\left({\mathit{\theta }}_D,{\mathit{\theta }}_A\right)$. 5:     According to Equation (23), obtain angle estimations ${\widehat{\mathit{\theta }}}_D^{(i+1)}$ and ${\widehat{\mathit{\theta }}}_A^{(i+1)}$. 6:     According to Equation (18), obtain the path gains ${\widehat{\rho }}^{\left(i+1\right)}$. 7:     Delete path l, when ${\widehat{\rho }}^{\left(i+1\right)}<{\rho }_{th}$. 8: until${∥{\rho }^{\left(i+1\right)}-{\rho }^{\left(i\right)}∥}_2^{}<{\epsilon }_{th}$, and $L^{\left(i\right)}=L^{\left(i+1\right)}.$ 9: $\widehat{\rho }={\widehat{\rho }}^{\left(end\right)}$, ${\widehat{\mathit{\theta }}}_D^{}={\widehat{\mathit{\theta }}}_D^{\left(end\right)}$, and ${\widehat{\mathit{\theta }}}_A^{}={\widehat{\mathit{\theta }}}_A^{\left(end\right)}$.

The estimations ${\widehat{\mathit{\theta }}}_D^{\left(i\right)}$ and ${\widehat{\mathit{\theta }}}_A^{\left(i\right)}$ become accurate until the new estimations are nearly equal to the previous estimations. The final estimations ${\widehat{\mathit{\theta }}}_D^{\left(i\right)}$, ${\widehat{\mathit{\theta }}}_A^{\left(i\right)}$, and ${\widehat{\rho }}^{\left(i+1\right)}$ of the jitter channel parameters converge to their true values with high probability. So, the super-resolution jitter channel estimation scheme can be realized by moving the estimations ${\widehat{\mathit{\theta }}}_D^{\left(i\right)}$ and ${\widehat{\mathit{\theta }}}_A^{\left(i\right)}$ from initial on-grid coarse estimations to their actual off-grid positions.

4. Numerical Results

Simulation results are provided under different conditions, verifying the performance of the proposed jitter channel estimation scheme. The simulation parameters in the jitter multi-UAV mmWave massive MIMO system with hybrid beamforming are $N_R=N_T=64$, $N_X=N_Y^{}=32$, $L=3$, $N_R^{RF}=N_T^{RF}=1$, $d=\lambda /2$, and $\lambda =0.5$. Assuming the path gains are Gaussian, ${\rho }_l\sim CN\left(0,{\sigma }_{\rho }^2\right)$, every element of $\mathbf{X}$ follows ${\mathrm{X}}_{\mathrm{i},\mathrm{j}}=\sqrt{\rho /N_{\mathrm{T}}}{e}^{j{\omega }_{i,j}}$, where $\rho$ is the transmitted power, ${\omega }_{i,j}\sim \left[0,2\pi \right)$. Then, $SNR=\rho {\sigma }_{\alpha }^2/{\sigma }_n^2$, and ${\sigma }_n^2$ is the noise variance. Compared with OMP-based [28] and AMP-based [29] channel estimation schemes, the performance of the proposed jitter channel estimation is shown in Figure 3 and Figure 4.

The NMSE is defined as the performance metric

$$\mathbf{NMSE}=\mathbb{E}\left[{∥\mathbf{H}-\widehat{\mathbf{H}}∥}^2/\phantom{{∥\mathbf{H}-\widehat{\mathbf{H}}∥}^2{∥\mathbf{H}∥}^2}{∥\mathbf{H}∥}^2\right].$$

(24)

We, respectively, compare the NMSE performance under the LOS and non-line-of-sight (NLoS) channels. In both cases, the NMSE performance based on the proposed scheme is better than that of the OMP-based and AMP-based channel estimation schemes. In the channel estimation schemes, respectively, based on the OMP and AMP algorithms, the signal reconstruction effect is reduced, and the error is increased by the presence of noise in Figure 3 and Figure 4. In comparison, in the proposed jitter channel estimation scheme, all jitter AoAs/AoDs are optimized, and the interference between different paths is eliminated. So, the NMSE is smaller than the OMP scheme at the same SNR. Whether UPA or ULA antenna type is used, the super-resolution channel estimation is realized by the proposed scheme.

Compared with the OMP-based and AMP-based schemes, the spectral efficiency of the proposed jitter scheme is given with $N_R=N_T=16$, as shown in Figure 5. The spectral efficiency is evaluated in the multi-UAV mmWave massive MIMO system with hybrid beamforming. The case with ideal CSI was adopted as the upper bound for performance comparison. It can be observed in Figure 5 that the proposed super-resolution channel estimation is able to approach the upper bound. This is because the jitter AoAs/AoDs resolutions of the proposed scheme are decided without codebook size or angle quantization. Therefore, it can be concluded that the proposed scheme can achieve super-resolution channel estimation.

We observe the jitter robust performance of the proposed scheme with $N_R=N_T=16$ in Figure 6. The $\mathbf{NMS}{\mathbf{E}}_R$-based robustness is defined as the performance metric:

$$\mathbf{NMS}{\mathbf{E}}_R=\mathbb{E}\left[{∥\mathbf{H}-{\widehat{\mathbf{H}}}_0∥}^2/\phantom{{∥\mathbf{H}-{\widehat{\mathbf{H}}}_0∥}^2{∥\mathbf{H}∥}^2}{∥\mathbf{H}∥}^2\right],$$

(25)

where ${\widehat{\mathbf{H}}}_0$ is the channel estimations without jitter obtained by the improved IR-based OMP and AMP methods, respectively. Compared with the OMP and AMP schemes, the proposed IR-based scheme is more robust in the jitter channel estimation under the $\mathbf{NMS}{\mathbf{E}}_R$ performance against SNR. Different from the commonly used performance indicator NMSE in Figure 3 and Figure 4, we observe the jitter-robust performance of the proposed scheme based on the $\mathbf{NMS}{\mathbf{E}}_R$-based robustness in Equation (25) to prove which parameters influence beam alignment issues. It is shown in Figure 6 that the proposed channel estimation scheme based on IR and the Padam algorithm can obtain the final estimations ${\widehat{\mathit{\theta }}}_D^{\left(i\right)}$, ${\widehat{\mathit{\theta }}}_A^{\left(i\right)}$, and ${\widehat{\rho }}^{\left(i+1\right)}$ of the jitter channel parameters, converging to their true values with high probability. Additionally, the improved scheme can be realized by moving the estimations ${\widehat{\mathit{\theta }}}_D^{\left(i\right)}$ and ${\widehat{\mathit{\theta }}}_A^{\left(i\right)}$ from the initial on-grid coarse estimations to their actual off-grid positions. Thus, the jitter channel estimation accuracy is significantly improved.

5. Conclusions

We proposed a robust super-resolution channel estimation scheme for the jitter multi-UAV system with hybrid beamforming. The jitter channel estimation problem is transformed as the optimization of objective function for sparsity-weighted summation and data-fitting error. Initializing the on-grid coarse points, the proposed scheme can iteratively find the off-grid actual positions by the Padam algorithm. Numerical results show that the improved scheme obviously increases the accuracy of the off-grid jitter AoAs/AoDs estimation. The robust performance of the proposed scheme is state-of-art, even in jitter channel estimation. In consideration of large-scale antennas for the multi-UAV mmWave massive MIMO system, more robust channel estimation models with reduced complexity could be studied in our future work.