1 Introduction

In recent years, there has been a rapid progress in the wireless communication [1–3], and many wireless transmission techniques have been proposed to meet the requirement of ultra-reliable and low-latency [4–7]. The unmanned aerial vehicle (UAV), an aerial platform developed with modern communication technologies, is an unmanned and reusable aircraft powered by electricity or fuel. The UAV has wide applications and can be used in intelligence investigation, disaster relief and rescue, precise target location and tracking, etc. The UAV can increase the coverage area of 5th generation (5G) network and increase the capacity of 5G systems and user equipments (UEs) as well. Therefore, it has received a lot of attention by the industry and academia.

UAV communications usually adopt the mmWave band, since the multi-input multi-output (MIMO) technology is easy to use in this band [8, 9]. The MIMO system, equipped with multiple antennas at both the transmitter and receiver, makes use of spatial resources adequately and has high data throughput and reliable communication quality. Channel tracking/prediction is always a research hot spot in MIMO systems. Channel state information (CSI) is indispensable for the data decoding at the receiver, which is often obtained by channel estimation. Usually, it is completed with the aide of the pilot signal sent by the transmitter frequently. However, by doing so, quite a part of the channel resources are occupied by the pilot signal, especially for fast fading channels. To reduce the length of the pilot signal on the whole, channel tracking or prediction was proposed by some literature [10]. These methods include the Kalman filter-based method [10], extended Kalman filter-based method [11], sequential Monte Carlo filter-based method [12], and particle filter-based method [13]. Recently, there are some works on the channel tracking in UAV MIMO systems. In [14], a channel tracking method for UAV MIMO communication systems was proposed and investigated, in which the method explores the characteristics of time-varying UAV channels with the beam squint effect. In [15], to improve the quality of the UAV navigation, the authors designed a channel tracking algorithm for its flight control system, where the time-varying spatial channel is characterized by a 3D geometry-based channel model and the algorithm incorporates the outputs of multiple sensors in order to reduce the training overhead and energy consumption.

On the other front, channel predictive information can be feeded back to the transmitter for preprocessing in order to improve the system performance in terms of the capacity or the bit-error-rate (BER) [16]. For instance, in [16], the authors proposed a subspace-based channel tracking scheme for precoded MIMO orthogonal frequency division multiplexing (MIMO-OFDM) systems. The predictive CSI was used for precoding at the transmitter to improve the system throughput. With the development of UAV communications, it can be foreknown that there will be a large amount of data to be transmitted interactively between UAVs and ground base stations/terminals. For instance, with the development of forestry informatization, forestry resource management needs more data support. The UAV is suitable for collecting forestry data and sending them back to the ground control center. Besides, the UAV can be used for data offloading of possibly overloaded cellular base stations in hot spots such as sports venues and cinemas. It can also be used for periodic data distribution/collection in large-scale Internet of Things networks. Large amount of data transmission requires temporary storage space. The caching technique pre-stores the data during the off-peak traffic and improves 5G system performance in many aspects [17–19]. With the aid of cache, it is convenient for the UAV to store data for its own use or forward data to other UAVs or ground terminals [20–22]. Clearly, the feedback information of channel prediction can also be used for UAV MIMO systems to improve the quality of data transmission.

Most previous work above is conducted based on the linear vector space, i.e., Euclidean space, which can be viewed as a linear manifold. Recently, with the rise of information geometry, a few channel tracking methods based on the non-linear manifold were proposed and investigated [23–26]. For instance, in [23], the authors first modeled the column spaces spanned by the right singular matrix of the MIMO channel as a point of Grassmannian manifold, and then proposed a method to track the movement of the point as well as an adaptive codebook for precoding. In [26], for the time-varying MIMO system, the authors proposed a predictive quantizer for the eigenvectors of the Gramian matrix that is created from the channel matrix, in which the quantizer operates on the compact Stiefel manifold. In this paper, we mainly predict a sub-matrix of the channel right singular matrix, which is created by selecting a few columns of the channel right singular matrix. The sub-matrix has more rows than columns and can be modeled as a point of Grassmannian manifold. When the channel is time varying, it is reasonable to assume that this sub-matrix moves along the geodesic within a short period of time. Therefore, we intend to use the theory and results of Grassmannian manifold in the channel prediction.

For the prediction methods based on Grassmannian manifold, a key parameter in the design is the choice of the search step. However, most previous work did not investigate the parameter thoroughly. Furthermore, considering the complexity of the practical UAV MIMO channel, this article incorporates both the time correlation and spatial correlation of the channel and studies the effects of various system parameters on the chordal error performance. The contributions of this article are briefly listed as follows:

  1. 1

    We propose a novel channel predictive precoding method which is based on the geodesic equation on the Grassmannian manifold and analyze its computational complexity as well as the conventional one. Compared with the conventional method, although the complexity of the two methods is the same, the proposed method has better prediction performance in terms of the channel chordal error.

  2. 2

    We adopt a general Ricean channel model which incorporates the time and spatial correlation by combining the Kronecker model and Jake’s model. The effects of various system parameters on the system performance in terms of the chordal error of channel predictor as well as the choice of optimal step are thoroughly investigated.

  3. 3

    In the Appendix, we prove a corollary which is inferred from the geodesic equation and is helpful to verify the correctness of the geodesic equation from two points on the Grassmannian manifold.

The rest of this paper is organized as follows. The closed-loop UAV-enabled MIMO system model is described in Section 2. The process of channel precoding prediction and the corresponding problem of interest are illustrated in Section 3. The effect of various system parameters on the performance of the prediction algorithm is presented in Section 4, followed by conclusions in Section 5.

Notations: Boldface lowercase letters denote vectors, and boldface uppercase letters denote matrices. The notation E(.) represents the statistical expectation; \(\mathbb {C}^{\textit {m}\times \textit {n}}\) and \(\mathbb {U}^{\textit {m}\times \textit {n}}\) represent the space consisting of m×n complex matrices and the space consisting of m×n orthogonal complex matrices, respectively. For a vector x, \( \mathbf {x} \sim {\mathbb {CN}} \left ({{\boldsymbol {\mu }},{\mathbf {R}}} \right) \) represents that x follows a complex Gaussian distribution with mean μ and covariance matrix R; for a matrix X, the notations X1/2, XH, and Tr(X) denote its square root, Hermitian transpose, and trace, respectively; besides, Im is an m×m identity matrix, Im,n with n<m is created by selecting the first n columns of Im, and Um is an m×m unitary matrix.

2 The closed-loop UAV-enabled MIMO system model

2.1 The UAV MIMO channel

As shown in Fig. 1, we consider a UAV-enabled MIMO system, where the base station (BS) and the UAV are equipped with NT and NR antennas, respectively [27, 28]. Both the BS and UAV are equipped with cache storage units. With the aid of cache, the UAV can store data conveniently for its own use or forward them to other UAVs/ground terminals at the next time. Figure 1a is the original system diagram, from which we see that there exists the line of sight (LOS) path between BS and UAV; due to the reflections of high buildings, trees, or others, there also exist non-LOS (NLOS) paths. Hence, the channel between BS and UAV is modeled as a Ricean fading channel.

Fig. 1
figure 1

MIMO system model with prediction feedback link

Full size image

For the UAV MIMO channel \({\mathbf {H}}(t)\in {\mathbb {C}}^{N_{R}\times N_T}\), we not only consider its time variation, but also the spatial correlation. In [29], a spatial-temporal correlated channel model was proposed; however, only transmitter correlation was considered. Herein, we also incorporate the receiver correlation of the channel. Consequently, the UAV MIMO channel is modeled as:

$$ {\mathbf{H}}\left(t \right) = {{{{\mathbf{H}}_{{\text{NLOS}}}}\left(t \right)} \over {\sqrt {{K_{Rice}} + 1} }}{\mathrm{ + }}\sqrt {{{{K_{Rice}}} \over {{K_{Rice}} + 1}}} {{\mathbf{H}}_{{\text{LOS}}}}\left(t \right). $$
(1)

In Eq. (1), the NLOS path HNLOS(t) is further expressed as:

$$ {{\mathbf{H}}_{{\text{NLOS}}}}\left(t \right) = {\boldsymbol{\theta }}_{R}^{1/2}{{\mathbf{H}}_{\omega} }\left(t \right){\boldsymbol{\theta }}_{T}^{1/2}, $$
(2)

where θT and θR denote the transmitter and receiver correlation matrix of the channel, respectively; KRice is the Ricean factor. \({{\mathbf {H}}_{\omega } }\left (t \right) \buildrel \Delta \over = \left ({{h_{ij}}\left (t \right)} \right) \in \mathbb {C}^{N_{R} \times N_T}\) and it is a random matrix that has i.i.d complex Gaussian random variables with zero mean and unit variance. The autocorrelation of each entry follows Jake’s model and it is given by:

$$ {\mathrm{E}}\left[ {{h_{ij}}\left({{t_1}} \right){{\left({{h_{ij}}\left({{t_2}} \right)} \right)}^*}} \right]{\mathrm{ = }}{J_0}\left({2\pi {f_d}\left({{t_2} - {t_1}} \right)} \right), $$
(3)

in which the function J0(.) is the well-known Bessel function of the first kind, fd is the Doppler shift, for the NLOS path, which is due to the motions of the UAV, scatters, etc.

Assuming that both BS and the UAV adopt the linear antenna array, the LOS path HLOS(t) is expressed as [30]:

$$ {{\mathbf{H}}_{{\text{LOS}}}}\left(t \right) = {e^{- j2\pi {{f'}_d}t}}{{\mathbf{x}}_R}{\left({{{\mathbf{x}}_T}} \right)^{H}}, $$
(4)
$$\begin{aligned} {{\mathbf{x}}_R}{\mathrm{ = }} \left[\begin{array}{cccc} 1&{{e^{- j2\pi \cos \left(\alpha \right){d_R}/\lambda }}}& \cdots &{{e^{- j2\pi \left({{N_R} - 1} \right)\cos \left(\alpha \right){d_R}/\lambda }}} \end{array}\right]^{T},\\ {{\mathbf{x}}_T}{\mathrm{ = }} \left[\begin{array}{cccc} 1&{{e^{- j2\pi \cos \left(\beta \right){d_T}/\lambda }}}& \cdots &{{e^{- j2\pi \left({{N_T} - 1} \right)\cos \left(\beta \right){d_T}/\lambda }}} \end{array}\right]^{T}, \end{aligned} $$

where fd′ is Doppler shift for the LOS path, which is due to the motion of the UAV; dR and dT are the antenna spacing at BS and the UAV, respectively; λ is the wavelength; and α and β are the angles of arrival and departure, respectively.

Note that the channel defined in Eq. (1) is assumed to be flat fading. If the paths are resolvable, we can incorporate the OFDM technique and transform the frequency selective channel into a flat fading one.

Remark: Concerning the angles of arrival and departure, the definitions are explained as follows. Take the angle of arrival as an example. As in Fig. 2, a plane electromagnetic wave arrives at the uniform linear antenna array. We adopt the definition of [30] and define α as the direction of arrival (DOA), i.e., the angle between the incident direction of the electromagnetic wave and the x-axis. Consequently, the time delay between the two adjacent antennas is dRcos(α)/λ, which determines the form of the steering vector xR. Whereas in [15], the authors define α′ as the DOA. Consequently, the time delay between the two adjacent antennas is dRsin(α′)/λ. Clearly, the two definitions are equivalent since sin(α′)=cos(α). Besides, it is worth noting that the latter requires the assumption that the incident direction of the wave is only in the x-y plane but the former does not require any assumptions. Therefore, we adopt the former definition.

Fig. 2
figure 2

The uniform linear antenna array and the arrival electromagnetic wave

Full size image

2.2 The system operation

Before transmission, the user signal \({\mathbf {x}}\left (t \right) \in \mathbb {C}^{N_{s}\times 1}\), sent by the transmitter, is multiplied by a precoding matrix \({\mathbf {W}}\left (t \right) \in {\mathbb {U}^{{N_T} \times {N_S}}}\), where t is the time index and NS is the number of data streams. After passing through fading channels, the signal at the receiver can be formulated as:

$$ {\mathbf{y}}\left(t \right) = {\mathbf{H}}\left(t \right){\mathbf{W}}\left(t \right){\mathbf{x}}\left(t \right) + {\mathbf{n}}\left(t \right), $$
(5)

in which x(t) has independent data streams with its autocorrelation function being \({\mathrm E}\left [ {{\mathbf {x}}\left (t \right){{\mathbf {x}}^{H}}\left (t \right)} \right ] = {{\mathbf {I}}_{{N_S}}}\), and \({\mathbf {n}}\left (t \right) \sim {\mathbb {CN}} \left ({\mathbf {0}},\sigma _{n}^{2} {\mathbf {I}}_{N_{R} \times N_R} \right)\) is the AWGN noise [31, 32], where the effect of noise on the communication systems can be found in [33, 34]. Notation \({\mathbf {H}}\left (t \right) \in \mathbb {C}^{N_{R} \times N_T} \) denotes the MIMO flat fading channel, where the fading scenarios can be found in [35–37].

Assume that the channel estimation is ideal and the receiver has perfect channel state information (CSI) [38–40], H(t). If the channel is quasi-static fading, which means that two adjacent channel matrices can be viewed as identical, i.e., H(t+1)≈H(t), the optimal precoder at time (t+1) is given by:

$$ {{\mathbf{W}}^{OPT}}\left({t + 1} \right) = {{\mathbf{V}}_d}\left(t \right){{\mathbf{U}}_d}, $$
(6)

where \({{\mathbf {U}}_d} \in \mathbb {U}^{{N_S} \times N_S}\) is an arbitrary unitary matrix and \({{\mathbf {V}}_d}\left (t \right) \in {\mathbb {U}^{{N_T} \times {N_S}}}\) is conducted by selecting the first NS columns of V(t), which comes from singular value decomposition (SVD) of H(t), expressed as:

$$ {\mathbf{H}}\left(t \right){\mathrm{ = }}{\mathbf{U}}\left(t \right){\mathbf{\Lambda }}\left(t \right){{\mathbf{V}}^{H}}\left(t \right). $$
(7)

The precoder expressed by (6) can not only maximize the system asymptotic capacity at high signal to noise ratio (SNR), but also minimize the mean square error of data detection.

However, if the channel is fast fading, the precoder at current time t is no longer optimal for the next time t+1. In this case, as in Fig. 1b, we need a predictor, which produces a predicted precoder \({{\mathbf {\widetilde {V}}}_d}\left (t + 1 \right)\) for time t+1 by the use of Vd(t) and Vd(t−1), where the definition of Vd(t−1) is similar to Vd(t). Consequently, the precoder at time t+1 is given by:

$$ {\mathbf{W}}\left({t + 1} \right) = {{\mathbf{\widetilde V}}_d}\left({t + 1} \right). $$
(8)

Clearly, the smaller the distance between \({{\mathbf {\widetilde V}}_d}\left ({t + 1} \right)\) and Vd(t+1), the better the prediction performance. We use the chordal distance to measure their distance, the details of which will be given in the next section.

In addition, in the process of channel prediction, there is a key parameter termed as the step parameter. An aim of this article is to study the effect of various channel settings on the optimal step parameter as well as the expected chordal distance (error).

3 Channel prediction precoding method

3.1 Preliminaries of the Grassmannian manifold

Grassmannian manifold, denoted as \({\mathbb {G}}\left ({{N_T},{N_S}} \right)\) with NS<NT, is important in mathematical differential geometry and has many engineering applications such as channel tracking [41, 42], image identification [43], and emotion recognition [44] in recent years. It is the set of all NS-dimensional linear subspaces of an n-dimensional space. The corresponding definition is given by:

$$ {\mathbb{G}}\left({{N_T},{N_S}} \right){\mathrm{ = }}\left\{ {\left[ {\mathbf{X}} \right] \in {\mathbb{C}^{{N_T} \times {N_S}}},{{\mathbf{X}}^{H}}{\mathbf{X}} = {{\mathbf{I}}_{{N_S}}}} \right\}, $$
(9)

where [X] is an equivalent class that is a set defined by \(\left [ {\mathbf {X}} \right ] \!\buildrel \Delta \over = \!\left \{\! {{\mathbf {X}}{{\mathbf {U}}_{{N_S}}}\!:\!{{\mathbf {U}}_{{N_S}}}\,{\text {is}}\,{\text {unitary}}}\! \right \}\). For the sake of convenience, we might as well write \({\mathbf {X}} \in {\mathbb {G}}\left ({{N_T},{N_S}} \right)\), meaning that X is an equivalent class whose columns span the same p-dimensional subspace [16, 45]. For two points \({{\mathbf {X}}_1},{{\mathbf {X}}_2} \in {\mathbb {G}}\left ({{N_T},{N_S}} \right)\), the chordal distance is usually used to measure their distance and it is defined as [14, 23]:

$$ D\left({{{\mathbf{X}}_1},{{\mathbf{X}}_2}} \right){\mathrm{ = }}{\sqrt {{N_S} - \left\| {{{\mathbf{X}}_1}^{H}{{\mathbf{X}}_2}} \right\|} _F}, $$
(10)

where ∥·∥F denotes the Frobenius norm.

Besides, the concept of geodesic plays an important role in information geometry. The term geodesic comes from geodesy, the science of measuring the Earth. Initially, a geodesic was the shortest route between two points on the Earth’s surface. Now, it is generalized and defined as the shortest route between two points of a manifold. The geodesic emanating from X1 to X2 is expressed as [41, 42]:

$$ {\mathbf{X}}\left(s \right) = {{\mathbf{Q}}_1} \exp \left({{\mathbf{B}}s} \right){{\mathbf{I}}_{{N_T},{N_S}}}, 0 \le s \le 1. $$
(11)

where \({{\mathbf {Q}}_1} \in {\mathbb {U}^{{N_T} \times {N_T}}} = \left ({{{\mathbf {X}}_1}\;{\mathbf {X}}_{1}^ \bot } \right)\), and \({\mathbf {X}}_{1}^ \bot \in {\mathbb {U}^{{N_T} \times \left ({{N_T} - {N_S}} \right)}}\) is some orthogonal basis of the orthogonal complement of X1. The matrix \({\mathbf {B}} \in {\mathbb {C}^{{N_T} \times {N_T}}}\) is skew Hermitian, and it has the following form:

$$ {\mathbf{B}}{\mathrm{ = }}\left({\begin{array}{*{20}{c}} {\mathbf{0}}&{ - {{\mathbf{A}}^{H}}}\\ {\mathbf{A}}&{\mathbf{0}} \end{array}}\right), $$
(12)

where \({\mathbf {A}} \in {\mathbb {C}^{\left ({{N_T} - {N_S}} \right) \times {N_S}}}\) is referred to as the velocity matrix. The matrix A can be further written as:

$$ {\mathbf{A}}{\mathrm{ = }}{{\mathbf{U}}_2}\Phi {\mathbf{U}}_{\mathrm{1}}^{H}, $$
(13)

where Φ is diagonal with its elements φi, 1≤i≤NS. The matrices \({\mathbf {U}}_{1} \in {\mathbb {U}^{{N_S} \times {N_S}}}\) and \({\mathbf {U}}_{2} \in {\mathbb {U}^{\left ({{N_T} - {N_S}} \right) \times {N_S}}}\) come from the cosine-sine decomposition of a specially constructed matrix [41, 42]:

$$ \left({\begin{array}{c} {{\mathbf{X}}_{1}^{H}{{\mathbf{X}}_2}}\\ {{{\left({{\mathbf{X}}_{1}^ \bot } \right)}^{H}}{{\mathbf{X}}_2}} \end{array}}\right) = \left({\begin{array}{cc} {{{\mathbf{U}}_1}}&{\mathbf{0}}\\ {\mathbf{0}}&{{{\mathbf{U}}_2}}\\ \end{array}}\right)\left({\begin{array}{c} {\mathbf{C}}\\ {\mathbf{S}} \end{array}}\right){\mathbf{V}}_{1}^{H}, $$
(14)

in which C is diagonal with its entries cosφi on the diagonal; S is also diagonal with its entries sinφi, on the diagonal, 1≤i≤NS; and V1 is the right singular matrix of the SVD of \({\mathbf {X}}_{1}^{H}{{\mathbf {X}}_2}\), which can also be inferred from Eq. (14).

With Eqs. (12)–(14), the key matrix B will be created and the geodesic Eq. 11 can be finally given afterwards. According to the original definition of X(s), we have the following corollary:

$$ {\mathbf{X}}\left({\mathrm{1}} \right) = {{\mathbf{Q}}_1}\exp \left({\mathbf{B}} \right){{\mathbf{I}}_{{N_T},{N_S}}}{\mathrm{ = }}{{\mathbf{X}}_2} $$
(15)

We believe that the above is important since it can verify the correctness of the geodesic equation; however, it was not proved in [41, 42]. Herein, we supply a proof which is shown in the Appendix.

3.2 Review of the conventional channel prediction method [41]

The matrix Vd(t) can be viewed as a point of Grassmannian manifold. It is reasonable to assume that the matrix Vd(t) moves or changes along the geodesic within a short period of time. The conventional prediction obtains the next point mainly by two steps: (1) parallel transport the tangent matrix at Vd(t−1) and (2) construct a new geodesic for prediction.

First, given two points Vd(t) and Vd(t−1), the tangent matrix at the point Vd(t−1), in the direction of the next point Vd(t), is given by [41]:

$$ {\mathbf{\Delta }}\left({t - 1} \right){\mathrm{ = }}{\mathbf{V}}_{d}^ \bot \left({t - {\mathrm{1}}} \right){\mathbf{A}}\left({t - 1} \right), $$
(16)

where \({\mathbf {V}}_{d}^ \bot \left ({t - {\mathrm {1}}} \right) \in {\mathbb {U}^{{N_T} \times \left ({{N_T} - {N_S}} \right)}}\) is some orthogonal basis of orthogonal complement of Vd(t−1), \({\mathbf {A}}\left ({t - 1} \right) \in {\mathbb {C}^{\left ({{N_T} - {N_S}} \right) \times {N_S}}}\) is the velocity matrix, and other details are easy to deduce by referring to the previous subsection.

Then, we transport the tangent matrix \({\mathbf {\Delta }}\left ({t - 1} \right) \in {\mathbb {C}^{{N_T} \times {N_S}}}\) parallel along the geodesic between Vd(t) and Vd(t−1), and obtain a new tangent matrix [45]:

$$ {\mathbf{\widetilde {\Delta}}}\left({t,{s_{\mathrm{0}}}} \right){\mathrm{ = }}\left[ \begin{array}{l} {{\mathbf{I}}_{{N_T}}} - {{\mathbf{U}}_F}{\mathbf{U}}_{F}^{H} + {{\mathbf{V}}_d}\left({t - 1} \right){{\mathbf{V}}_F}\sin \left({{{\mathbf{\Lambda }}_F}{s_{\mathrm{0}}}} \right){\mathbf{U}}_{F}^{H}\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + {{\mathbf{U}}_F}\cos \left({{{\mathbf{\Lambda }}_F}{s_{\mathrm{0}}}} \right){\mathbf{U}}_{F}^H \end{array}\right]{\mathbf{\Delta }}\left({t - 1} \right), $$
(17)

where 0≤s0≤1 is the step parameter and \({{\mathbf {U}}_F}{{\mathbf {\Lambda }}_F}{\mathbf {V}}_{F}^{H}\) is the compact SVD of the matrix Δ(t−1). It is clear that \({\mathbf {\widetilde {\Delta }}}\left ({t,{\mathrm {0}}} \right) = {\mathbf {\Delta }}\left ({t - 1} \right)\) and \({\mathbf {\widetilde \Delta }}\left ({t,{\mathrm {1}}} \right)\) is the tangent matrix at the point Vd; we shall transport the tangent matrix Δ(t−1) parallel from Vd(t−1) to Vd(t), which means s0=1.

Finally, the predictive point on the manifold at time t+1, \({\mathbf {\overline {V}}}_{d}^{{\text {Conv}}}\left ({t,s} \right)\), can be derived as [23, 41, 45]:

$$ {\mathbf{\overline {V}}}_{d}^{{\text{Conv}}}\left({t,s} \right) = {{\mathbf{V}}_d}\left(t \right){{\mathbf{V}}_E}\cos \left({{{\mathbf{\Lambda }}_E}s} \right){\mathbf{V}}_{E}^{H} + {{\mathbf{U}}_E}\sin \left({{{\mathbf{\Lambda }}_E}s} \right){\mathbf{V}}_{E}^{H}, $$
(18)

where \({{\mathbf {U}}_E}{{\mathbf {\Lambda }}_E}{\mathbf {V}}_{E}^{H}\) is the compact SVD of the matrix \({\mathbf {\widetilde {\Delta } }}\left ({t,{\mathrm {1}}} \right)\) and s is the step parameter which can be further optimized. The optimal step sOPT can be found by solving the following equation:

$$ {s_{OPT}} = \mathop {\arg \min } \limits_{s} {\mathrm E} \left[ {D\left({{\mathbf{\overline {V}}}_{d}^{{\text{Conv}}}\left({t,s} \right),{{\mathbf{V}}_d}\left({t + 1} \right)} \right)} \right]. $$
(19)

Once the optimal step is found, we substitute it into (18) to obtain \({\mathbf {\widetilde {V}}}_{d}^{{\text {Conv}}}\left ({t + 1} \right) \buildrel \Delta \over = {\mathbf {\overline {V}}}_{d}^{{\text {Conv}}}\left ({t,{s_{OPT}}} \right)\).

3.3 The proposed channel prediction method

In this subsection, we propose a novel channel prediction method, in which the next point is predicted by use of the geodesic from Vd(t−1) to Vd(t) directly. First, according to the previous subsection A, the geodesic emanating from Vd(t−1) to Vd(t) is expressed as:

$$ {{\mathbf{F}}}\left({t - 1,s} \right) = {{\mathbf{Q}}_1}\left({t - 1} \right)\exp \left[ {{\mathbf{B}}\left({t - 1} \right)s} \right]{{\mathbf{I}}_{{N_T},{N_s}}} $$
(20)

where \({{\mathbf {Q}}_1}\left ({t - 1} \right) \in {\mathbb {U}^{{N_T} \times {N_T}}} = \left ({{\mathbf {V}}_{d} \left ({t - {\mathrm {1}}} \right) {\mathbf {V}}_{d}^ \bot \left ({t - {\mathrm {1}}} \right)} \right)\); \({\mathbf {V}}_{d}^ \bot \left ({t - {\mathrm {1}}} \right)\) and B(t−1) have similar definitions in the previous subsection A. Clearly, F(t−1,0)=Vd(t−1) and F(t−1,1)=Vd(t).

Second, the predictive point on the manifold at time t+1 is given by:

$$ {\mathbf{F}}\left({t - 1,s + 1} \right)\; = {{\mathbf{Q}}_1}\left({t - 1} \right)\exp \left[ {\left({s + 1} \right){\mathbf{B}}\left({t - 1} \right)} \right]{{\mathbf{I}}_{{N_T},{N_s}}},\;0 \le s \le 1. $$
(21)

Similarly, the optimal step sOPT in the above formula can be found by solving Eq. (19). Note that the function \({\mathbf {\overline {V}}}_{d}^{{\text {Conv}}}\left ({t,s} \right)\) in Eq. (19) shall be replaced with F(t−1,s+1). Once it is found, the precoding matrix or predictive precoder is solved by \({\mathbf {\widetilde {V}}}_{d}^{{\text {Prop}}}\left ({t + 1} \right) \buildrel \Delta \over = {\mathbf {F}}\left ({t - 1,{s_{OPT}} + 1} \right)\).

The complete prediction process is summarized by Algorithm 1.

Note that in the calculation of the matrix function exp[B(t−1)s], we shall use the following matrix series [46]:

$$ \exp \left[ {{\mathbf{B}}\left({t - 1} \right)s} \right]{\mathrm{ = }}\sum\limits_{i} {{{\left[ {{\mathbf{B}}\left({t - 1} \right)s} \right]}^{i}}{\mathrm{/}}i!} \,. $$
(22)

3.4 Computational complexity

Usually, the number of floating point flops is used to measure of the complexity of an algorithm. We define a flop as one floating point operation and it has computational complexity O(1). Note that the matrix product X1X2 requires O(mnk) flops, where \({\mathbf X}_{1} \in \mathbb C^{m\times n}\) and \({\mathbf X}_{2} \in \mathbb C^{n\times k}\); for an n×n matrix, both its inverse and eigen-decomposition operations require O(n3) flops; for an m×n matrix with m≥n, the complexity of its singular value decomposition (SVD) is O(m2n) if the complete left singular matrix is needed; otherwise, it is O(mn2) [47].

To begin with, we analyze the computational complexity of the proposed method. Initially, the two points Vd(t−1) and Vd(t) are generated from the SVD of H(t−1) and H(t), respectively. Since the left singular matrices are not needed, this operation has complexity of \(O\left (n_{max}n_{min}^{2}\right)\), where nmax=max(NT,NR),nmin=min(NT,NR). Then, it is found that step 1 requires \(O\left (N_{T}^{2} N_{S}\right)\) flops and step 2 requires \(O\left (N_{T}^{3}\right)\) flops; step 3 requires \(O\left (N_{S}^{2} N_{T}\right) + O(N_{S}^3) +O\left [\left (N_T-N_{S}\right)N_{T} N_S)\right ]=O\left (N_{T}^{2} N_{S}\right)\) flops, since NS≤NT; and step 4 requires \(O\left (N_T-N_S)N_{S}^{2}\right)\) flops. For the calculation of exp[B(t−1)s] in step 5, by using the techniques in the Appendix, we have that:

$$\exp \left[ {{\mathbf{B}}\left({t \,-\, 1} \right)s} \right] \,=\, \left({\begin{array}{ccc} {{{\mathbf{U}}_1}\left({t - 1} \right)\cos \left({\Phi}s \right){\mathbf{U}}_{1}^{H}(t-1)}\! &\!{ - {{\mathbf{U}}_1}\left({t - 1} \right)\sin \left({\Phi}s \right){\mathbf{U}}_{2}^{H}(t-1)}\\ {{{\mathbf{U}}_2}\left({t - 1} \right)\sin \left({\Phi}s \right){\mathbf{U}}_{1}^{H}\left({t - 1} \right)} \!&\!{{{\mathbf{U}}_2}\left({t - 1} \right)\cos \left({\Phi}s \right){\mathbf{U}}_{2}^{H}\left({t - 1} \right)} \end{array}}\right). $$

With the above, given s, solving F(t−1,s) has complexity of \(O\left (N_{T} N_{S}^{2} + N_{T}^{2} N_{S}\right)=O\left (N_{T}^{2} N_{S}\right)\). Assuming that the minimal search interval that stands for the accuracy of sOPT is ξ, 1/ξ times of calculating F(t−1,s) are required to find sOPT. Hence, step 5 has complexity of \( O\left (N_{T}^{2} N_S/ \xi \right)\). In summary, the complexity of the proposed method is \(O\left (n_{max}n_{min}^{2}\right)+ O\left (N_{T}^{3}\right)+O\left (N_{T}^{2} N_S/ \xi \right)=O\left [\left (n_{max}n_{min}^{2}\right)+ N_{T}^{2} \left (N_T+N_S/ \xi \right)\right ]\), where nmax=max(NT,NR),nmin=min(NT,NR).

Then, we analyze the complexity of the conventional method briefly. Initially, we also need the two points Vd(t−1) and Vd(t). In step 1, to generate Vd⊥(t−1), assume that the same method as in the subsection C is used. Further, we note that for an m×n matrix with m≥n, the compact SVD has complexity of O(mn2) because the full right singular matrix is not needed. Paying attention to the above points, it is easy to find that the complexity of the conventional method is the same as the proposed method.

4 Results and discussion

In this section, computer simulation is deployed to investigate the performance of the prediction algorithm. The UAV MIMO channel samples are generated according to the channel model in Section 2. Specifically, a 4×4 UAV MIMO system with NS = 2 is adopted. For the transmitter and receiver correlation matrices, θT and θR, the well-known exponential correlation model is adopted [48], in which the (i,j)th entry of the correlation matrix is ρ|i−j|, and the constant ρ is the spatial correlation coefficient. We further denote the correlated coefficients at the transmitter and receiver as ρT and ρR, respectively. The angles of arrival and departure α and β are set to 30∘ and 20∘, respectively; the antenna spacing dT=dR=λ. Besides, the Doppler shifts for the LOS path and NLOS path are set to be identical, i.e., fd′=fd.

4.1 Performance comparison between the proposed method and the conventional one

Figure 3 compares the chordal error performance of the conventional method with that of the proposed one, where ρT = 0.2 and ρR = 0.3. Two cases are considered: fdT=0.02 and KRice=0 dB for case 1 and fdT=0.05 and KRice=10 dB for case 2, where fdT is the normalized Doppler shift, fd is the Doppler shift, and T is the sample interval. Observe that given s>0, the chordal error for the proposed method is always less than that for the conventional method. For case 1 in Fig. 3a, the minimum chordal error for the conventional method is about −15.42 dB at s=0.9; the minimum chordal error for the proposed one is −18.01 dB at s=1.0 and it is less than −15.42 dB. Therefore, the proposed method is superior to the conventional method. Similar phenomena can be observed in Fig. 3b.

Fig. 3
figure 3

Performance comparison between the conventional method and the proposed one

Full size image

Figure 4 compares two methods with different values of normalized Doppler shift fdT, in which ρT = 0.2, ρR = 0.3, and KRice = 0 dB. Observe that for both methods, given the step parameter s, increasing fdT results in increasing chordal error. Increasing fdT also has the optimal step parameter sOPT decrease. For instance, for the proposed method, the optimal steps are about 1, 0.9, and 0.6, for fdT =0.01, 0.05, and 0.1, respectively. Similar to Fig. 3, the proposed method has better performance than the conventional method, no matter whether the normalized Doppler shift is large or not.

Fig. 4
figure 4

Two methods’ performance comparison with different values of normalized Doppler shift fdT

Full size image

4.2 The effects of a few key system parameters

Figure 5 shows the optimal step parameter sOPT with different normalized Doppler shifts. The system settings are consistent with those in Fig. 4. Observe that for increasing normalized Doppler shift, sOPT decreases from 1 to 0.6 gradually. The optimal step sOPT is found by solving Eq. (19), where the searching step is set to 0.1, and hence, the curve shows a stepwise downward trend.

Fig. 5
figure 5

The optimal step sOPT with different values of normalized Doppler shift fdT

Full size image

Figures 6 and 7 present the effects of ρT and ρR on the chordal error, respectively. For both figures, fdT = 0.05 and KRice = 0 dB. Further, we set ρR = 0.2 for Fig. 6 and ρT = 0.2 for Fig. 5, respectively. From Fig. 6, we observe that given the step parameter s, increasing ρT will decrease the chordal error and also decrease the optimal s. For instance, given s = 0.5, the chordal errors for ρT = 0, 0.2, 0.4, 0.6, and 0.8 are −9.2 dB, −9.8 dB, −11.0 dB, −12.8 dB, and −15.8 dB, respectively; the optimal s decreases from 0.9 to 0.8 slightly when increasing ρT from 0 to 0.8. In Fig. 7, increasing the correlation coefficient also reduces the chordal error, but conversely, the optimal s increases from 0.8 to 0.9 slightly with increasing ρR. Therefore, these results indicate that both ρT and ρR are negatively related to the chordal error; ρT is slightly negatively related to the optimal step parameter, whereas the opposite is true for ρR.

Fig. 6
figure 6

The effect of transmitter correlated coefficient ρT on the chordal error

Full size image
Fig. 7
figure 7

The effect of receiver correlated coefficient ρR on the chordal error

Full size image

Figure 8 depicts the effect of the Ricean factor KRice on the chordal error, where fdT = 0.05, ρT = 0.2, and ρR = 0.3. Observe that given s, larger KRice results in smaller chordal error. This is possibly because that the increase of KRice reduces the randomness of the UAV channel variation and tracking the LOS path is easier than the NLOS path. In addition, the optimal step sOPT is about 0.9 and nearly unaffected by KRice.

Fig. 8
figure 8

The effect of the Ricean factor KRice on the chordal error

Full size image

5 Conclusions

This article studies the problem of channel predictive precoding in UAV MIMO systems. A novel predictive precoding method that is based on the geodesic in the Grassmannia manifold is proposed. Instead of transporting the tangent vector parallel along the geodesic, the proposed method predicts the next point by use of the geodesic equation directly. We compare the proposed method with the conventional one by simulation and analyze their computational complexity. Besides, the effects of various system parameters, including time and spatial correlated coefficients and the Ricean factor, on the chordal error of the predictor are also thoroughly investigated.

The results reveal the following. First, the normalized Doppler shift is positively related to the chordal error, whereas the transmitter and receiver correlation coefficients and the Ricean factor are negatively related to the chordal error. Second, the normalized Doppler shift is negatively related to the optimal step parameter; the transmitter correlation coefficient is slightly negatively related to the optimal step parameter, whereas the opposite is true for the receiver correlation coefficient. Third, with the same computational complexity, the proposed predictive algorithm is superior to the conventional method in terms of the channel predictive error—the channel chordal error. The above research results will provide a reference for possible system implementation in the future.

A further study is to investigate the effect of channel estimation errors on the predictor, where the two channel estimates include the previous and current estimates which are not assumed to be ideal. How to design a robust predictor is worthy of researching. Another possible study is to incorporate the codebook quantizer into the predictor design.

6 Appendix

6.1 Proof of Eq. (15)

First, the matrix function exp(B) can be expanded as:

$$\begin{aligned} \exp \left({\mathbf{B}} \right)&{\mathrm=}\sum\limits_{i} {{{\mathbf{B}}^{i}}{\mathrm{/}}i!} \\ &= {{\mathbf{I}}_N} + \left({\begin{array}{cc} {\mathbf{0}}&{ - {{\mathbf{U}}_1}\Phi {\mathbf{U}}_{\mathrm{2}}^{H}}\\ {{{\mathbf{U}}_2}\Phi {\mathbf{U}}_{\mathrm{1}}^{H}}&{\mathbf{0}} \end{array}}\right)/1!\; + \left({\begin{array}{cc} { - {{\mathbf{U}}_1}{\Phi ^{2}}{\mathbf{U}}_{1}^{H}}&{\mathbf{0}}\\ {\mathbf{0}}&{ - {{\mathbf{U}}_2}\Phi {^{2}}{\mathbf{U}}_{2}^{H}} \end{array}}\right)/2!\\ &+ \left({\begin{array}{cc} {\mathbf{0}}&{{{\mathbf{U}}_1}{\Phi ^{3}}{\mathbf{U}}_{2}^{H}}\\ { - {{\mathbf{U}}_2}{\Phi ^{3}}{\mathbf{U}}_{1}^{H}}&{\mathbf{0}} \end{array}}\right)/3!\; + \left({\begin{array}{cc} {{{\mathbf{U}}_1}{\Phi ^{4}}{\mathbf{U}}_{1}^{H}}&{\mathbf{0}}\\ {\mathbf{0}}&{{{\mathbf{U}}_2}{\Phi ^{4}}{\mathbf{U}}_{2}^{H}} \end{array}}\right)/4!\\ &+ \left({\begin{array}{cc} {\mathbf{0}}&{ - {{\mathbf{U}}_1}{\Phi ^{5}}{\mathbf{U}}_{2}^{H}}\\ {{{\mathbf{U}}_2}{\Phi ^{5}}{\mathbf{U}}_{1}^{H}}&{\mathbf{0}} \end{array}}\right)/5!\; + \;\; \cdots ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~(A-1) \end{aligned} $$

Since

$$ \sin \Phi {\mathrm{ = }}{\sum\nolimits}_{i = 1}^{+ \infty } {{{\left({ - 1} \right)}^{i}}} {\Phi^{2i{\mathrm{ + }}1}}/\left({2i + 1} \right)! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{(A-2)} $$

and

$$ \cos \Phi {\mathrm{ = }}{\sum\nolimits}_{i = 0}^{+ \infty } {{{\left({ - 1} \right)}^{i}}} {\Phi^{2i}}/\left({2i} \right)!, ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{(A-3)} $$

the formula (A−1) can be simplified as:

$$ \exp \left({\mathbf{B}} \right) = \left(\!\! {\begin{array}{ccc} {\mathbf{0}}&{ - {{\mathbf{U}}_1}\sin \Phi {\mathbf{U}}_{2}^{H}}\\ {{{\mathbf{U}}_2}\sin \Phi {\mathbf{U}}_{1}^{H}}&{\mathbf{0}} \end{array}}\!\!\right) \!+ \! \left(\!\! {\begin{array}{cc} {{{\mathbf{U}}_1}\cos \Phi {\mathbf{U}}_{1}^{H}}&{\mathbf{0}}\\ {\mathbf{0}}&{{{\mathbf{U}}_2}\cos \Phi {\mathbf{U}}_{2}^{H}} \end{array}}\!\!\right). ~~~~~~~~~~~~~~{(A-4)} $$

Then, we have:

$$ \begin{aligned} {\mathbf{X}}\left({\mathrm{1}} \right) &= {{\mathbf{Q}}_1}\exp \left({\mathbf{B}} \right){{\mathbf{I}}_{{N_T},{N_S}}}\\ &{\mathrm{ = }}{{\mathbf{Q}}_1}\left\{ {\left(\!\! {\begin{array}{cc} {{{\mathbf{U}}_1}\cos \Phi {\mathbf{U}}_{1}^{H}}&{\mathbf{0}}\\ {\mathbf{0}}&{{{\mathbf{U}}_2}\cos \Phi {\mathbf{U}}_{2}^{H}} \end{array}}\!\!\right)} \right.\left. {{\mathrm{ \!\,+\,\! }}\left(\!\!{\begin{array}{cc} {\mathbf{0}}&{ - {{\mathbf{U}}_1}\sin \Phi {\mathbf{U}}_{2}^{H}}\\ {{{\mathbf{U}}_2}\sin \Phi {\mathbf{U}}_{1}^{H}}&{\mathbf{0}} \end{array}}\!\!\right)} \right\}{{\mathbf{I}}_{{N_T},{N_S}}}\\ &{\mathrm{ = }}\left({{{\mathbf{X}}_1}\;\!{\mathbf{X}}_{1}^ \bot } \right) \!\! \left\{\!\! {\left(\!\! {\begin{array}{cc} {{{\mathbf{U}}_1}\cos \Phi {\mathbf{U}}_{1}^{H}}\!\!\!\!&\!\!\!\!{\mathbf{0}}\\ {\mathbf{0}}&{{{\mathbf{U}}_2}\cos \Phi {\mathbf{U}}_{2}^{H}} \end{array}}\!\! \right)} \right.\left. {{\mathrm{ \!\!\!+ \!\! }}\left(\!\! {\begin{array}{cc} {\mathbf{0}}\!\!\!&\!\!\!{ - {{\mathbf{U}}_1}\sin \Phi {\mathbf{U}}_{2}^{H}}\\ {{{\mathbf{U}}_2}\sin \Phi {\mathbf{U}}_{1}^{H}}\!\!\!\!&\!\!\!\!{\mathbf{0}} \end{array}}\!\!\right)} \!\!\right\}\!{{\mathbf{I}}_{{N_T},{N_S}}}\\ &{\mathrm{ = }}\left[ {\left({\begin{array}{ccc} {{{\mathbf{X}}_1}{{\mathbf{U}}_1}\!\cos \Phi {\mathbf{U}}_{1}^{H}}&{{\mathbf{X}}_{1}^ \bot {{\mathbf{U}}_2}\!\cos \Phi {\mathbf{U}}_{2}^{H}} \end{array}}\right)} \right.\left. {{\mathrm{+}} \left({\begin{array}{cc} {{\mathbf{X}}_{1}^ \bot {{\mathbf{U}}_2}\!\sin \Phi {\mathbf{U}}_{1}^{H}}\!\!&\!\!{ - {{\mathbf{X}}_1}\!{{\mathbf{U}}_1}\!\sin \Phi {\mathbf{U}}_{2}^{H}} \end{array}}\right)} \right]\!{{\mathbf{I}}_{{N_T},{N_S}}}\\ &= {{\mathbf{X}}_1}{{\mathbf{U}}_1}\cos \Phi {\mathbf{U}}_{1}^{H} + {\mathbf{X}}_{1}^ \bot {{\mathbf{U}}_2}\sin \Phi {\mathbf{U}}_{1}^H \end{aligned} ~~~~~~{(A-5)} $$

Further manipulation of (A−5) results in:

$$ \begin{aligned} &{\mathbf{X}}\left({\mathrm{1}} \right) = \left({{{\mathbf{X}}_1}\;{\mathbf{X}}_{1}^ \bot } \right)\left({\begin{array}{*{20}{c}} {{{\mathbf{U}}_1}}&{\mathbf{0}}\\ {\mathbf{0}}&{{{\mathbf{U}}_2}} \end{array}}\right)\left({\begin{array}{cc} {\cos \Phi }\\ {\sin \Phi } \end{array}}\right){\mathbf{U}}_{1}^{H}\\ &= \left({{{\mathbf{X}}_1}\;{\mathbf{X}}_{1}^ \bot } \right)\left({\begin{array}{ccc} {{{\mathbf{U}}_1}}&{\mathbf{0}}\\ {\mathbf{0}}&{{{\mathbf{U}}_2}} \end{array}}\right)\left({\begin{array}{cc} {\mathbf{C}}\\ {\mathbf{S}} \end{array}}\right){\mathbf{U}}_{1}^{H}\\ &= \left({{{\mathbf{X}}_1}\;{\mathbf{X}}_{1}^ \bot } \right)\left({\begin{array}{cc} {{{\mathbf{U}}_1}}&{\mathbf{0}}\\ {\mathbf{0}}&{{{\mathbf{U}}_2}} \end{array}}\right)\left({\begin{array}{cc} {\mathbf{C}}\\ {\mathbf{S}} \end{array}}\right){\mathbf{V}}_{1}^{H}{{\mathbf{V}}_1}{\mathbf{U}}_{1}^H \end{aligned} ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{(A-6)} $$

With (14), we have:

$$ \begin{aligned} {\mathbf{X}}\left({\mathrm{1}} \right){\mathrm{ = }}\left({{{\mathbf{X}}_1}\;{\mathbf{X}}_{1}^ \bot } \right)\left({\begin{array}{*{20}{c}} {{\mathbf{X}}_{1}^{H}{{\mathbf{X}}_2}}\\ {{{\left({{\mathbf{X}}_{1}^ \bot } \right)}^{H}}{{\mathbf{X}}_2}} \end{array}}\right){{\mathbf{V}}_1}{\mathbf{U}}_{1}^{H}\\ \;\;\;\;\;\;\;\;\;{\mathrm{ = }}\left({{{\mathbf{X}}_1}\;{\mathbf{X}}_{1}^ \bot } \right)\left({\begin{array}{*{20}{c}} {{\mathbf{X}}_{1}^{H}}\\ {{{\left({{\mathbf{X}}_{1}^ \bot } \right)}^{H}}} \end{array}}\right){{\mathbf{X}}_2}{{\mathbf{V}}_1}{\mathbf{U}}_{1}^{H}{\mathrm{ = }}{{\mathbf{X}}_2}{{\mathbf{V}}_1}{\mathbf{U}}_{1}^H \end{aligned} ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{(A-7)} $$

According to the definition of Grassmannia manifold, \({{\mathbf {X}}_2}{{\mathbf {V}}_1}{\mathbf {U}}_{1}^{H}\) is equivalent to X2. Therefore, Eq. (15) holds.