Tensor-Train Decomposition Based Hybrid Beamforming for Millimeter-Wave Massive MIMO/FSO in UAV with RIS Networks

1. Introduction

UAV (Unmanned aerial vehicle) combined with cellular networks can support communications in a cost-effective and highly mobile manner. UAV communication offers the possibility to establish new dedicated terrestrial networks [1-3]. In areas with dense traffic hotspots in cities, the coverage of traditional cellular base stations cannot meet the demand. UAV can quickly provide coverage and service triage to hotspots based on their high flexibility and easy deployment and have a wide range of application scenarios. Compared with terrestrial base stations, the ABS (Aerial Base Station) is more adaptable to environmental changes [4]. Therefore, UAV can also be deployed to provide emergency communication connections in areas without infrastructure coverage. Ground base stations are often destroyed during natural disasters such as earthquakes, tsunamis, and flash floods. The communication infrastructure in the disaster area is damaged and cannot provide communication services. This implements rescue operations more hindered.

To improve the transmission performance of UAV communication, RA-FSO-RF (Relay-Assisted Free Space Optical-Radio Frequency) was introduced relaying into the communication system [5-6]. The LDPC-MIMO/FSO (Low-Density Parity-Check Coded Multiple-Input Multiple-Output/Free Space Optical) method investigates the performance evaluation of power-series based málaga distributed MIMO/FSO links with m-QAM and pointing error (PE) [7]. DHPL-ML (Dynamic Hybrid Precoding with Manifold Learning) method proposes a manifold learning inspired dynamic hybrid precoding dual hop hybrid FSO-RF system based on antenna partitioning algorithm [8]. Similarly, MIMO techniques were used to eliminate the effects of atmospheric turbulence (AT) and pointing errors on link performance in FSO systems [9]. The bit error rate (BER) and the channel capacity of MIMO system employing an optical space scheme ware investigated under the combined effect of AT and PE [10]. The communication environment in UAV was greatly affected by atmospheric turbulence. It can have a significant impact on the performance of laser communication. The UAV communication mainly includes visual axis errors and jitter errors. At present, the jitter influence of free-space laser links was studied [11]. However, the existing algorithms have rarely reported the bit error rate performance of UAV under the influence of turbulence. Meantime, these previous studies on FSO-RF systems have focused on the FSO link or relay selection problem [12-15]. The problems of the PE and beam misalignment of millimeter-wave massive MIMO/FSO in UAV communication RF links are ignored.

The coverage performance of UAV base stations is affected by the dynamic movement of UAV and obstacles [16]. In addition, the issues of beam misalignment and Doppler frequency shift caused by obstacle interference and wind speed is needed to consider [17]. To improve UAV communication performance, the beamforming of UAV antenna arrays can help overcome the problems of beam misalignment and Doppler shift. Beamforming was designed for UAV communication with random blockages [18-19]. The 3D location of the UAV base station with a certain number of indoor users and obstacles was determined. UAV can perceive user and obstacle locations to maximize the total coverage of base stations. The ZF-HBF (Zero Forcing-Hybrid beamforming) method was explored for multi-user orthogonal frequency division multiplexing systems over downlink [20-23]. The zero forcing method is used to achieve beamforming according to the system's SE (Spectrum Effectiveness) up to the maximum upper bound. To enable fast beamforming training and tracking, a hierarchical structure of beamforming codebooks and design of hierarchical codebooks with different beam widths via sub-array techniques was investigated [24]. The Doppler effect as a result of UAV movement was examined and found that the Doppler effect may be catastrophic when low gain beamforming is used. A hybrid beamforming for millimeter-wave massive MIMO with channel estimation was proposed [25-27]. By exploiting the angular sparsity of mmWave channels, the continuously distributed angle of arrivals/departures (AoAs/AoDs) can be jointly estimated for hybrid beamforming. A Hy-BD (Hybrid Block Diagonalization) method on phase beamforming with a combiner was obtained [28]. The digital block diagonalization of massive array gain was used to process equivalent baseband channels. However, in the case of high-speed motion, the currently available beamforming techniques are limited in terms of performance improvement. The RIS (Reconfigurable Intelligent Surface) for beamforming in wireless communication [29] was considered for UAV communication systems. The RIS is needed to further consider the altitude, transmission distance, and multipath channel of the UAV to plan its path. The REHS (RIS for nonlinear energy harvesting algorithm) method was investigated to improve the UAV communication performance by optimizing the phase of the RIS reflective elements and the parameters of the nonlinear energy harvesting circuit [30-32]. However, the indoor wireless network coverage lacks adaptivity due to the mobility of UAV and the time-varying channel. The RBOP (RIS-Based optimization and passive beamforming algorithm) was investigated to maximize the signal received power of the UAV communication system [33]. The RBOP aims at joint trajectory design and passive beamforming in UAV communication. The system performance will degrade if the hybrid beamforming is targeted for partial connections. Multi-user beamforming for millimeter-wave massive MIMO in UAV leads to high-dimensional operations[34-37]. To maximize the spectral efficiency of the system, it is necessary to jointly design its RIS reflection coefficient and hybrid beamforming. These algorithms require high computational costs, leading to degraded energy efficiency.

To overcome the light PE and beam misalignment of millimeter-wave massive MIMO/FSO caused by UAV jitter, a Tensor-train decomposition based hybrid beamforming for FSO and millimeter-wave massive MIMO of UAV with RIS is investigated. Firstly, the RIS-assisted millimeter-wave massive MIMO/FSO in UAV channel model is established, and the high-dimensional channels of the RIS-assisted millimeter-wave massive MIMO/FSO of UAV is represented as the low-dimensional channels by Tensor-train decomposition. Secondly, the problem of maximizing the system spectral efficiency is decomposed into two sub-problems optimizing the light PE of FSO and the hybrid beamforming matrix. The FSO PE caused by UAV jitter can be effectively solved by BIGRU (Bidirectional Gated Recurrent Unit)-attention neural network model. The fast-fading channels and Doppler shift are estimated by the FCTPM (Fast Circulant Tensor Power Method) [38-39] based on the Tensor-train decomposition. Finally, the RIS phase shift matrix is optimized by the SVD (Singular Value Decomposition). The Hybrid beamforming and RIS phase shift matrix are estimated by the low-complexity PE-AltMin (Phase Extraction Alternating Minimization) method to solve the beam misalignment. Simulation experiments demonstrate that the proposed method has higher spectrum utilization than other methods.

2. System model

This paper considers a RIS-assisted multi-user millimeter-wave massive MIMO /FSO UAV communication system, as shown in Figure 1. Covering densely indoors can affect the connectivity of signal transmission, and the window of RIS components is small. When the UAV experiences shaking, it can cause beam misalignment during signal transmission. The multi-user millimeter-wave massive MIMO in UAV, in turn, leads to high dimensional operations in the beamforming process. To maximize the spectral efficiency of the system, it is necessary to jointly design its RIS reflection coefficient and hybrid beamforming. However, due to the non-convexity of the problem, it is decomposed into two independent optimization subproblems. The block diagram of a mixed FSO-RF system is presented in Figure 1 in which the source node (S) communicates with the destination node (D) through a decode-and-forward (DF) relay node (R). The S-R link is equipped with a single antenna, and the R-D link is a millimeter wave massive MIMO system model. Node R has both optical and RF signal processing capabilities. We employ non-coherent intensity modulation with direct detection (IM/DD) receiver at R. After converting the incoming optical signal to an electrical signal, node R utilizes a power splitter to separate the alternating current (AC) and direct current (DC) components. The unsolicited DC component (which is normally filtered out at the receiver) is applied to the energy harvesting unit which supplies the harvested power to the RF transmitter. The information-bearing AC component is given to the decoder circuit which decodes the information and demodulates using an RF modulation scheme before forwarding it through the RF transmitter.

2.1. FSO system model

The FSO channel between S and R is modelled as Gamma-Gamma distribution with pointing error [40]. The statistical behavior of the received optical irradiance is characterized by means of Gamma-Gamma turbulence model which has been widely used to model the FSO channel in the recent literature depending on its doubly stochastic scintillation model. The probability distribution function (PDF) of channel coefficient $h_{FSO}$ is given as

$$f_{h_{FSO}}\left(h_{FSO}\right)=\frac{{\varsigma }^2}{h_{FSO}\mathsf{\Gamma }\left(\alpha \right)\mathsf{\Gamma }\left(\beta \right)}{G}_{\mathrm{1,3}}^{\mathrm{3,0}}\left(\alpha \beta \frac{h_{FSO}}{I_{FSO}{M}_o}{|}_{{\varsigma }^2,\alpha ,\beta }^{{\varsigma }^2+1}\right)$$

(1)

Where $\mathsf{\Gamma }($.) is the well-known Gamma function,$G_{p,q}^{m,n}\left({\left.x\right|}_{b_1,\dots ,b_q}^{a_1,\dots ,a_p}\right)$is the Meijer’s G-function. The atmospheric turbulence induced fading channel gains of FSO links, denoted as $I_{FSO}$ deled with s and path loss. It is assumed that the elements of $I$ are modeled as independent and identically distributed random variables $\varsigma ={\overline{\mathsf{\omega }}}_e/2{\sigma }_s$ is the ratio of equivalent beam radius and zero boresight displacement standard deviation at the photo-detector (PD), the constant term $M_o$ is the power fraction that the detector receives when there is no, and $1/\mathsf{\alpha }$ and $1/\mathsf{\beta }$ are the variances of the large and small turbulence eddies, respectively. The FSO system is modeled as shown in Figure 2.

The parameters $\alpha$ and $\beta$ are the distance-dependent fading variables which corresponds to atmospheric turbulence conditions as

$$\alpha ={\left\{\mathrm{e}\mathrm{x}\mathrm{p}\left[\frac{0.49{\sigma }_R^2}{{\left(1+1.11{\sigma }_R^{\frac{12}{5}}\right)}^{\frac{7}{6}}}\right]-1\right\}}^{-1}$$

(2)

and

$$\beta ={\left\{\mathrm{e}\mathrm{x}\mathrm{p}\left[\frac{0.51{\sigma }_R^2}{{\left(1+0.69{\sigma }_R^{\frac{12}{5}}\right)}^{\frac{5}{6}}}\right]-1\right\}}^{-1}$$

(3)

Where ${\sigma }_R^2=1.23C_n^2{v}^{7/6}{d}_{SR}^{11/6}$ denoting the Rytov variance, $C_n^2$, $v$ and $d_{SR}$ represent the refractive index structure constant, wave number, and the link length, respectively. $C_n^2$ usually takes values in the range ${10}^{-17}-{10}^{-13}$ for weak up to strong turbulence conditions.

2.2. DF transmission protocol

The communication between S and D is accomplished in two time slots $T_1$ and $T_2$ for the first and second hop, respectively. The relay harvests energy during $T_1$ only as simultaneous harvesting and discharging of power increases the complexity at the node. The S employs sub-carrier intensity modulation (SIM) to convert RF signal vector $r$ with electrical power $\rho$ to an optical signal. A DC bias $A\in [A_{\mathrm{m}\mathrm{i}\mathrm{n}},A_{max}]$ (where $A_{\mathrm{m}\mathrm{i}\mathrm{n}}$ and $A_{max}$ are minimum and maximum values of DC bias, respectively) is added to the RF signal to ensure a non-negative optical signal. Let $P_{\mathrm{s}}$ represents the electrical power of S for transmitting the optical signal vector ${\mathrm{s}}_1$, then we can write:

$$s_1=\sqrt{P_{\mathrm{S}}}[\delta r+A]$$

(4)

Where $\delta$ is electrical to optical conversion coefficient. To prevent clipping due to non-linearity of laser diode such that it operates in linear region, $\delta$ should satisfy the following constraint [26]

$$\delta \le min\left(\frac{A-A_{min}}{\rho },\frac{A_{max}-A}{\rho }\right)$$

(5)

The electrical signal at the output of the PD, can be expressed as

$$y_{FSO}=h_{FSO}s+n_{FSO}$$

(6)

Where $h_{FSO}=\left(P_{\mathrm{S}}^2{\eta }^2/{\sigma }_{FSO}^2\right)I_{FSO}$ is the channel coefficient of S-D link, where $\eta$ is optical to electrical conversion coefficient. The noise ${\eta }_{FSO}$ is due to circuit noise as well as high-intensity background illumination and is conventionally modeled as being zero mean additive white Gaussian noise. ${\sigma }_{FSO}^2$is the variance of additive white Gaussian noise with zero mean.

2.3. Mathematical modeling of the PE tracking system

The received electrical SNR with the slow fading channel can be expressed by the following equation:

$$\mathrm{S}\mathrm{N}\mathrm{R}\left(h_{FSO}\right)=\frac{2P_t^2{h}^2{R}^2}{{\sigma }_n^2}$$

(7)

Where $h_{FSO}=h_p{h}_d$, it comprises the impact of PE $h_p$ which is an independent random variable, and $h_d$，the path loss component due to other atmospheric attenuation. The path loss component can be determined by the Beer-lambert's law as [33]

$$h_d=\mathrm{e}\mathrm{x}\mathrm{p}(-{\alpha }_{\mathrm{a}\mathrm{t}\mathrm{m}\mathrm{o}}z)$$

(8)

Where $z$ is the propagation distance and ${\alpha }_{\mathrm{a}\mathrm{t}\mathrm{m}\mathrm{o}}$ is the attenuation coefficient due to atmospheric conditions. The received power collected by the detector in the presence of geometric spread and PE can be expressed by the following equation:

$$h_p\approx A_0\mathrm{e}\mathrm{x}\mathrm{p}\left(\frac{-2r^2}{w_{\mathrm{z}\mathrm{e}\mathrm{q}}^2}\right)$$

(9)

Where $A_0=[\mathrm{e}\mathrm{r}\mathrm{f}(v){]}^2$ is the power collected at radial displacement of $\mathrm{r}=0$ and $w_{\mathrm{z}\mathrm{e}\mathrm{q}}^2\widehat{=}{w}_z^2\left(\right(\sqrt{\pi }\mathrm{e}\mathrm{r}\mathrm{f}\left(v\right))/(2v\mathrm{e}\mathrm{x}\mathrm{p}(-v^2)$is the equivalent beamwidth. $v=\sqrt{\pi }a/\sqrt{2}{w}_z$ is the beamwidth at distance $z$, $a$ is the radius of the detection aperture. Performance of FSO communication under PE can be evaluated based on the outage probability.

The performance of the proposed FSO system is provided in terms of outage probability due to s. The outage 0ccurs when the channel capacity, $C\left(h\right)$ is not sufficient to support the data rate $R_0$ and it can be expressed by the following equation:

$$P_{\mathrm{o}\mathrm{u}\mathrm{t}}\left(R_0\right)=P_r\left(C\right({\mathrm{S}\mathrm{N}\mathrm{R}}_{Rx}\left(h\right))<R_0)$$

$$P_{\mathrm{o}\mathrm{u}\mathrm{t}}\left(R_0\right)=P_r\left({\mathrm{S}\mathrm{N}\mathrm{R}}_{Rx}\left(h\right)<C^{-1}\left(R_0\right)\right)$$

(10)

$$=P_r(\frac{2P_t^2{h}^2{R}^2}{{\sigma }_n^2}<C^{-1}(R_0\left)\right).$$

(11)

The outage probability can be expressed by the following equation:

$$P_{\mathrm{o}\mathrm{u}\mathrm{t}}=P_r(h<h_0)$$

(12)

Where $h_0=\left(\right(\left(C^{-1}\right(R_0\left){\sigma }_n^2\right)/\left(2P_t^2{R}^2\right)){)}^{1/2}$. The outage probability is the cumulative density function of $h$ evaluated $h_0$ and it can be expressed by the following equation:

$$\begin{gathered} P_{\mathrm{o}\mathrm{u}\mathrm{t}}={\int }_0^{h_0}{f}_{h_p}\left(h_p\right)d{h}_p \\ =k^{\mathrm{\text{'}}}{\int }_0^{h_0}{(\mathrm{l}\mathrm{n}(\frac{A_0}{h_p}\left)\right)}^{m-1}(h_p{)}^{\zeta -1}d{h}_p. \end{gathered}$$

(13)

On substitution $x=\mathrm{l}\mathrm{n}(A_0/h_p)$

$$P_{\mathrm{o}\mathrm{u}\mathrm{t}}=k^{\mathrm{\text{'}}}{A}_0^{\zeta }{\int }_{\mathrm{\infty }}^{\ln \left(\frac{A_0}{h_0}\right)}{x}^{m-1}\mathrm{e}\mathrm{x}\mathrm{p}(-x\zeta )dx$$

(14)

Where

$$k^{\mathrm{\text{'}}}=\frac{{\zeta }^m}{\mathsf{\Gamma }\left(m\right)A_0^{{\zeta }^m}}$$

(15)

Using the approximation of

$${\int }_u^{\mathrm{\infty }}{x}^{v-1}\mathrm{e}\mathrm{x}\mathrm{p}(-\mu x)dx={\mu }^{-v}\mathsf{\Gamma }(v,u\mu )$$

(16)

$$P_{\mathrm{o}\mathrm{u}\mathrm{t}}=k^{\mathrm{\text{'}}}{A}_0^{\zeta }{\zeta }^{-m}\mathsf{\Gamma }(m,\zeta \mathrm{l}\mathrm{n}(\frac{A_0}{h_0}\left)\right)$$

(17)

2.4. Estimation of FSO PE with BIGRU-Attention model

When the distance of airborne laser communication is long, the requirements for the pointing of the airborne laser communication terminal are higher. Large s cause large deviations from the laser beam. The application of BiGRU-Attention neural network model in FSO system can effectively solve the problem of light PE caused by UAV jitter. The structure of BiGRU-Attention model is shown in Figure 3. The $y_{\mathrm{F}\mathrm{S}\mathrm{O}}$ in Equation (6), PE in Equation (9) and outage probability in Equation (17) are inputted into the neural network as parameters corresponding to $x_1$,$x_2$ and$x_3$, respectively. After the BiGRU-Attention layer, the BiGRU network extracts the time-varying information of the $y_{\mathrm{F}\mathrm{S}\mathrm{O}}$, PE and outage probability parameters. The attention mechanism will assign different attention weights to $y_{\mathrm{F}\mathrm{S}\mathrm{O}}$, PE and outage probability parameters. Different weights can distinguish the importance level of different parameters and improve the accuracy of classification. Finally, the classification result is obtained, so that the problems in FSO can be solved effectively.

The BiGRU network is a bi-directional recurrent neural network that simultaneously reads the inputs from both ends of the $y_{\mathrm{F}\mathrm{S}\mathrm{O}}$, PE, and outrage probability parameter time series This allows the BiGRU network to capture the temporal information of the above three parameter sequences of the inputs, which helps in extracting the feature information of the entire parameter. The attention mechanism network selects the most critical part of the $y_{\mathrm{F}\mathrm{S}\mathrm{O}}$, PE, outrage probability parameter time series by learning the attention distribution. This distribution indicates that the network should enhance the critical parameter feature information while decoding. The BiGRU-Attention model proposed in this paper enters the network through the input and output signals $y_{\mathrm{F}\mathrm{S}\mathrm{O}}$, PE and outage probability as parameters. Firstly, the BiGRU network learns the parameter time series information. Secondly, the importance of key parameter features is emphasized by the Attention network. Finally, the BiGRU-Attention network effectively solves the pointing error problem in the FSO system.

2.5. Multi-user UAV communication channel model with RIS

Firstly, the total number of users is K, and the number of data streams is $N_{s,K}$. The UAV communication transmitter is equipped with $N_T$ two transmitting antennas with $N_T^{RF}$ RF chains, where the $N_{s,K}\le N_T^{RF}\le N_T$. The RIS then consists of ${\mathcal{T}}^{\mathrm{‘}}$ reflective elements, where the set of reflective elements ${\mathcal{T}}^{\mathrm{‘}}$ =$\left\{1,\cdots ,{\mathrm{T}}^{\mathrm{‘}}\right\}$. In general, the RIS can communicate with the UAV using a specific link to exchange information about the channel status of the system, while allowing for better coordination of downlink transmissions to multipleuser sides. The reflection coefficient of the RIS is ${\varrho }_{t^{\mathrm{‘}}}=e^{j{\theta }_{t^{\mathrm{‘}}}},\forall t^{\mathrm{‘}}ϵ{\mathcal{T}}^{\mathrm{‘}}$.${ \theta }_{t^{\mathrm{‘}}}ϵ[\mathrm{0,2}\pi ）$denotes the first $t^{\mathrm{‘}}$ phase shift of the first RIS reflection element. The RIS reflection matrix is expressed as $\mathbf{\Phi }=diag\left\{{\varrho }_{1^{\mathrm{‘}}},\cdots ,{\varrho }_{{\mathcal{T}}^{\mathrm{‘}}}\right\}$, and that $\left|{\varrho }_{t^{\mathrm{‘}}}\right|=1,\forall t^{\mathrm{‘}}ϵ{\mathcal{T}}^{\mathrm{‘}}$. Similarly, the phase shifts of the RIS reflection elements have discrete values satisfying ${\theta }_{t^{\mathrm{‘}}}ϵ\mathcal{I}\triangleq \left\{\frac{2\pi \stackrel{⃑}{i}}{2^{\stackrel{⃑}{I}}}|\stackrel{⃑}{i}=\mathrm{0,1},\cdots ,2^{\stackrel{⃑}{I}}-1\right\}$, and $\mathcal{I}$ denotes the resolution of each reflecting element.

Let $s_1=\sqrt{P_{\mathrm{R}}}{y}_{FSO}$ denotes the signal forwarded by the node R. For this system, for K terrestrial users, the signal is first passed through a baseband beamforming matrix ${\mathit{F}}_{BB}$ into a digital signal vector. Then the simulated beamforming matrix ${\mathit{F}}_{\mathrm{R}\mathrm{F}}$ gets the UAV transmit signal as:

$$\mathit{x}={\sum }_{k=1}^K{\mathit{F}}_{RF}{\mathit{F}}_{BB,k}{{\mathit{s}}_1}_k$$

(18)

Where ${\mathit{F}}_{RF}ϵ{\mathbb{C}}^{N_T\times N_{RF}^T}$,${\mathit{F}}_{BB,k}=\left[{\mathit{F}}_{BB,1},\cdots {\mathit{F}}_{BB,K}\right]ϵ{\mathbb{C}}^{N_{RF}^T\times N_{s,k}}$ ,$N_{s,k}={\sum }_{k=1}^K{N}_{s,k}$.${{\mathit{s}}_1}_k={\left[{\mathit{s}}_1^T,\cdots ,{\mathit{s}}_K^T\right]}^Tϵ{\mathbb{C}}^{N_s\times 1}$ denotes$k$ the concatenation of the parallel data streams of the individual users, and ${\mathit{s}}_kϵ{\mathbb{C}}^{N_{S,k}\times 1}$ denotes the first $k$ Gaussian encoded data stream vector for the first user satisfying $\mathbb{E}\left[{\mathit{s}}_k{\mathit{s}}_k^H\right]={\mathit{I}}_{N_{s,k}}$. The maximum signal transmission power $P$ satisfies ${‖{\mathit{F}}_{RF}{\mathit{F}}_{BB}‖}_F^2\le P$.

The UAV transmits a communication signal to the RIS and the RIS passes the communication signal to theuser. Theuser is equipped with $N_R$ a receiving antenna with $N_R^{RF}$ an RF chain, where ${N_{S,k}\le N}_R^{RF}\le N_R$ . For signal accuracy, an analogue combiner matrix is designed ${\mathit{W}}_{RF,k}$ and a digital combiner matrix ${\mathit{W}}_{BB,k}（{\mathit{W}}_{RF,k}ϵ{\mathbb{C}}^{N_R\times N_{RF}^R}$ and ${\mathit{W}}_{BB,k}ϵ{\mathbb{C}}^{N_{RF}^R\times N_{s,k}}$ ). Therefore, the received signal ${\mathit{y}}_k\in {\mathbb{C}}^{{\mathsf{{\rm N}}}_{s,k}\times 1}$ at the kth terrestrial user can be expressed by the following equation:

$${\mathit{y}}_k\left(t\right)={\mathit{W}}_{BB,k}^H{\mathit{W}}_{RF,k}^H（{\ddot{\mathit{M}}}_k\left(t\right){\mathbf{\Phi }\ddot{\mathit{R}}\left(t\right)）\mathit{F}}_{RF}{\mathit{F}}_{BB,k}{\mathit{s}}_k$$

$$+{\mathit{W}}_{BB,k}^H{\mathit{W}}_{RF,k}^H{\left({\ddot{\mathit{M}}}_k\left(t\right)\mathbf{\Phi }\ddot{\mathit{R}}\left(t\right)\right)\sum _{j\ne k,jϵ\mathcal{K}}{\mathit{W}}_{BB,k}^H{\mathit{W}}_{RF,k}^H{{\mathit{F}}_{RF}\mathit{F}}_{BB,k}{\mathit{s}}_k+{\mathit{W}}_{BB,k}^H{\mathit{W}}_{RF,k}^H\mathit{n}}_k\left(t\right)$$

(19)

Where ${\mathit{W}}_{BB,k}^H{\mathit{W}}_{RF,k}^H（{\ddot{\mathit{M}}}_k\left(t\right){\mathbf{\Phi }\ddot{\mathit{R}}\left(t\right)）\mathit{F}}_{RF}{\mathit{F}}_{BB,k}{\mathit{s}}_k$ represents the received signal excluding other interference,${\mathit{W}}_{BB,k}^H{\mathit{W}}_{RF,k}^H{\left({\ddot{\mathit{M}}}_k\left(t\right)\mathbf{\Phi }\ddot{\mathit{R}}\left(t\right)\right)\sum _{j\ne k,jϵ\mathcal{K}}{\mathit{W}}_{BB,k}^H{\mathit{W}}_{RF,k}^H{{\mathit{F}}_{RF}\mathit{F}}_{BB,k}{\mathit{s}}_k}_{}$is the interference signal and ${\mathit{W}}_{BB,k}^H{\mathit{W}}_{RF,k}^H{\mathit{n}}_k$ represents the total noise.${\mathit{F}}_{BB,k}ϵ{\mathbb{C}}^{N_{RF}^T\times N_{s,k}}$ and ${\mathit{F}}_{RF}ϵ{\mathbb{C}}^{N_T\times N_{RF}^T}$ represent the baseband beamforming matrix and the analogue beamforming matrix, respectively. ${\mathit{n}}_k~\mathcal{C}\mathcal{N}\left(0,{\delta }^{2}I\right)$ is the $k$ additive Gaussian white noise at theuser. ${\ddot{\mathit{M}}}_k\left(t\right)$ is the millimeter wave large-scale channel matrix from the RIS to theuser, and $\ddot{\mathit{R}}\left(t\right)$ is the millimeter-wave large-scale channel matrix from the UAV to the RIS. ${\ddot{\mathit{M}}}_k\left(t\right)$ and $\ddot{\mathit{R}}\left(t\right)$ are defined as follows, respectively:

$${\ddot{\mathit{M}}}_k\left(t\right)=\sqrt{\frac{{{\mathrm{T}}^{\mathrm{‘}}N}_R}{Q_{{\ddot{M}}_k}}}\sum _{q=1}^{Q_{{\ddot{M}}_k}}{\stackrel{̿}{\kappa }}_{km,q}\left(t\right){\stackrel{̿}{\mathit{\alpha }}}_R\left({{\phi }_{km,q}}^R\left(t\right)\right){\stackrel{̿}{\mathit{\alpha }}}_{P_1}^H\left({{\phi }_{km,q}}^T\left(t\right),{{\eta }_{km,q}}^T\left(t\right)\right)ϵ{\mathbb{C}}^{N_R\times {\mathrm{T}}^{\mathrm{‘}}}$$

(20)

$$\ddot{\mathit{R}}\left(t\right)=\sqrt{\frac{{{\mathrm{T}}^{\mathrm{‘}}N}_T}{Q_{\ddot{R}}}}\sum _{q=1}^{Q_{\ddot{R}}}{\stackrel{̿}{\kappa }}_{r,q}\left(t\right){\stackrel{̿}{\mathit{\alpha }}}_{P_2}\left({{\phi }_{r,q}}^R\left(t\right),{{\eta }_{r,q}}^R\left(t\right)\right){\stackrel{̿}{\mathit{\alpha }}}_T^H\left({{\nu }_{r,q}}^T\left(t\right)\right)ϵ{\mathbb{C}}^{{\mathrm{T}}^{\mathrm{‘}}\times N_T}$$

(21)

The UAV and theuser are equipped with uniform linear arrays. RIS is equipped with a uniform planar array of size ${{\mathrm{T}}^{\mathrm{‘}}}_y\times {{\mathrm{T}}^{\mathrm{‘}}}_z$. $Q_{{\ddot{M}}_k}$ and $Q_{\ddot{R}}$ are the numbers of channel paths from the RIS to the $k$user and the number of channel paths between the UAV and the RIS. ${\stackrel{̿}{\kappa }}_{km,q}\left(t\right)$ denotes the path gain from the RIS to the $k$user , the path gain between ${\stackrel{̿}{\kappa }}_{r,q}\left(t\right)$ denotes the fast fading path gain between the UAV and the RIS.${\stackrel{̿}{\mathit{\alpha }}}_R（{{\phi }_{km,q}}^R\left(t\right))$,${\stackrel{̿}{\mathit{\alpha }}}_T（{{\nu }_{r,q}}^T\left(t\right)$, ${\stackrel{̿}{\mathit{\alpha }}}_{P_1}（{{\phi }_{km,q}}^T\left(t\right)$,${{\eta }_{km,q}}^T\left(t\right)$ and ${\stackrel{̿}{\mathit{\alpha }}}_{P_2}（{{\phi }_{r,q}}^R\left(t\right),{{\eta }_{r,q}}^R\left(t\right)）$are the steering vectors corresponding to the uniform linear array and the uniform planar array, which can be defined, respectively:

$${\stackrel{̿}{\mathit{\alpha }}}_R\left({{\phi }_{km,q}}^R\left(t\right)\right)=\frac{1}{\sqrt{N_R}}\left[1,\dots ,e^{j\left(N_R-1\right)d\sin \left({{\phi }_{km,q}}^R\left(t\right)\right)\frac{2\pi }{\lambda }}\right]$$

(22)

$${\stackrel{̿}{\mathit{\alpha }}}_T\left({{\nu }_{r,q}}^T\left(t\right)\right)=\frac{1}{\sqrt{N_T}}\left[1,\dots ,e^{j\left(N_T-1\right)d\sin \left({{\phi }_{km,q}}^T\left(t\right)\right)\frac{2\pi }{\lambda }}\right]$$

(23)

$${\stackrel{̿}{\mathit{\alpha }}}_{P_1}\left({{\phi }_{km,q}}^T\left(t\right),{{\eta }_{km,q}}^T\left(t\right)\right)=\frac{1}{\sqrt{M_y{M}_z}}\left[1,\dots ,e^{jd\sin \left(\left(M_y-1\right)\sin \left({{\phi }_{km,q}}^T\left(t\right)\right)\sin \left({{\eta }_{km,q}}^T\left(t\right)\right)+\left(M_z-1\right)\cos \left({{\phi }_{km,q}}^T\left(t\right)\right)\right)\frac{2\pi }{\lambda }}\right]$$

(24)

$${\stackrel{̿}{\mathit{\alpha }}}_{P_2}\left({{\phi }_{r,q}}^R\left(t\right),{{\eta }_{r,q}}^R\left(t\right)\right)=\frac{1}{\sqrt{M_y{M}_z}}\left[1,\dots ,e^{jd\sin \left(\left(M_y-1\right)\sin \left({{\phi }_{r,q}}^R\left(t\right)\right)\sin \left({{\eta }_{r,q}}^R\left(t\right)\right)+\left(M_z-1\right)\cos \left({{\eta }_{r,q}}^R\left(t\right)\right)\right)\frac{2\pi }{\lambda }}\right]$$

(25)

where $d=\frac{\lambda }{2}$ is the antenna spacing and$\lambda$ is the signal wavelength.

3. Power fading factor and Doppler shift estimation for fast fading channels based on Tensor-train Decomposition

3.1. Tensor-train Decomposition

Tensor-train is a higher-order tensor decomposition, which represents a higher-order tensor as a product of several lower-order tensor. Tensor-train decomposition is an effective method for compressed representation of high-dimensional data, which can reduce the number of dimensions and simplify the computation, while maintaining the main structure of the data. It is very useful in dealing with large-scale high-dimensional data. The main benefit of the Tensor-train method is the ability to obtain a low-dimensional representation of the high-dimensional signal by Tensor-train decomposition, which simplifies the dimensionality and improves the robustness. Thus Tensor-train helps to improve the accuracy of parameter estimation for power fading and Doppler shift.

$\mathcal{X}$ is an nth order tensor of size ${\mathit{I}}_1\times \cdots \times {\mathit{I}}_N$. Its Tensor-train representation uses N tensor cores ${\mathcal{G}}_1,{\mathcal{G}}_2,\cdots ,{\mathcal{G}}_N$ where the size of each tensor core is ${\mathit{R}}_{n-1}\times {\mathit{I}}_n\times {\mathit{R}}_n$, ${\mathit{R}}_0={\mathit{R}}_N=1$.

$$\mathcal{X}=\sum _{r_1=1}^{R_1}\sum _{r_2=2}^{R_2}\cdots \sum _{r_{N-1}=1}^{R_{N-1}}{\mathcal{G}}_1\left(:,r_1\right)\circ {\mathcal{G}}_2\left(r_1,:,r_2\right)\circ \cdots \circ {\mathcal{G}}_{\mathrm{N}-2}\left(r_{N-2},:,r_{\mathrm{N}-1}\right)\circ {\mathcal{G}}_{\mathrm{N}}\left(r_{N-1},:\right)$$

(26)

A Tensor-train $\mathcal{Y}={\mathcal{G}}_1\circ {\mathcal{G}}_2·\cdots \circ {\mathcal{G}}_N$ of rank $\left(R_0,R_1,\cdots R_N\right)$ has rank $R_1{R}_2\cdots R_{N-1}$, and the equivalent Kruskal tensor is denoted as $\mathcal{Y}=⟦{\mathit{A}}^{\left(1\right)},{\mathit{A}}^{\left(2\right)},\cdots {\mathit{A}}^{\left(N\right)}⟧$, where ${\mathit{A}}^{\left(n\right)}$ is given by the following equation:

$${\mathit{A}}^{\left(n\right)}={\left[{\mathcal{G}}_n\right]}_{\left(2\right)}\left({\mathit{I}}_{R>n}^T⨂{\mathit{I}}_{R_{n-1}{R}_n}⨂{\mathit{I}}_{R<n-1}^T\right)$$

(27)

Where $R_{<n}=R_0{R}_1\cdots R_{n-1}$， $R_{>n}=R_{n+1}\cdots R_{N-1}{R}_{N }$.

3.2. Estimation of Power Fading Factor and Doppler Shift for Fast Fading Channels using Tensor-Train Decomposition

The received signal at the single-user of a millimeter-wave MIMO in UAV can be expressed as the following equation:

$${\mathbf{y}}_{\mathit{k}}\left(t\right)={\mathit{W}}_{BB,k}^H{\mathit{W}}_{RF,k}^H\sqrt{\frac{{{\mathrm{T}}^{\mathrm{‘}}N}_R}{Q_{{\ddot{M}}_k}}}\sum _{q=1}^{Q_{{\ddot{M}}_k}}{\stackrel{̿}{\kappa }}_{km,q}\left(t\right){\stackrel{̿}{\mathit{\alpha }}}_R\left({{\phi }_{km,q}}^R\left(t\right)\right){\stackrel{̿}{\mathit{\alpha }}}_{P_1}^H\left({{\phi }_{km,q}}^T\left(t\right),{{\eta }_{km,q}}^T\left(t\right)\right){\mathit{F}}_{RF}{\mathit{F}}_{BB,k}{\mathit{s}}_k$$

$$+{\left({\mathit{W}}_{RF}{\mathit{W}}_{BB}\right)}^H\overleftarrow{\mathit{N}}\left(t\right)$$

(28)

Set the conversion again as follows:$\sqrt{\frac{N_T{N}_R}{Q}}$ is a constant, namely $\sqrt{\frac{N_T{N}_R}{Q}}=1$.

$${\tilde{\mathit{\alpha }}}_R\left({{\phi }_{km,q}}^R\left(t\right)\right)={\mathit{W}}_{BB,k}^H{\mathit{W}}_{RF,k}^H{\stackrel{̿}{\mathit{\alpha }}}_R\left({{\phi }_{km,q}}^R\left(t\right)\right)ϵ{\mathbb{C}}^{N_{RF}^R\times 1}$$

(29)

$${\tilde{\mathit{\alpha }}}_P\left({{\eta }_{km,q}}^T\left(t\right)\right)={\stackrel{̿}{\mathit{\alpha }}}_{P_1}^H\left({{\phi }_{km,q}}^T\left(t\right),{{\eta }_{km,q}}^T\left(t\right)\right){\ddot{\mathbf{\Phi }\mathit{R}}\left(t\right)\mathit{F}}_{RF}{\mathit{F}}_{BB,k}{\mathit{s}}_kϵ{\mathbb{C}}^{N_{RF}^T\times 1}$$

(30)

Where ${\tilde{\mathit{\alpha }}}_R\left({{\phi }_{km,q}}^R\left(t\right)\right)$ and ${\tilde{\mathit{\alpha }}}_P\left({{\eta }_{km,q}}^T\left(t\right)\right)$ represent the receive vector and the send vector respectively. Equation (28) is expressed as the following equation:

$$\mathbf{Y}\left(t\right)=\sum _{q=1}^{Q_{{\ddot{M}}_k}}{\stackrel{̿}{\kappa }}_{km,q}\left(t\right){\tilde{\mathit{\alpha }}}_R\left({{\phi }_{km,q}}^R\left(t\right)\right){\tilde{\mathit{\alpha }}}_P\left({{\eta }_{km,q}}^T\left(t\right)\right)+{\left({\mathit{W}}_{RF}{\mathit{W}}_{BB}\right)}^H\overleftarrow{\mathit{N}}\left(t\right)$$

(31)

The received signal tensor for T consecutive time slots as a rank one vector accumulation pattern is expressed as the following equation:

$$\mathcal{Y}=\sum _{q=1}^Q\mathcal{T}{\tilde{\mathit{\alpha }}}_R\left({{\beta }_q}^R\left(t\right)\right)·{\tilde{\mathit{\alpha }}}_T\left({{\gamma }_q}^T\left(t\right)\right)·{\tilde{\kappa }}_q\left(t\right)+\mathcal{N}$$

(32)

Where $\mathcal{T}$ is the weight of the multiplication of the rank one factors,$\mathcal{Y}\mathsf{ϵ}{\mathbb{C}}^{N_{RF}^R\times N_{RF}^T\times T}$ is the received signal tensor, ${\tilde{\kappa }}_q\left(t\right)={\left[{\kappa }_q\left(1\right),\cdots ,{\kappa }_q\left(T\right)\right]}^Tϵ{\mathbb{C}}^{T\times 1}$ is the vector containing the fast fading coefficients from T consecutive time slots, and $\mathcal{N}ϵ{\mathbb{C}}^{N_{RF}^R\times N_{RF}^T\times T}$ is the noise tensor. Meanwhile, equations (32) can continue to be expressed as a Tensor-train decomposition of the$K$ $ruskal$form:

$$\mathcal{Y}=\mathcal{T};{\mathit{U}}_1,{\mathit{U}}_2,{\mathit{U}}_3+\mathcal{N}$$

(33)

$${\mathit{U}}_1=\left[{\tilde{\mathit{\alpha }}}_R\left({{\mathit{\beta }}_1}^R\left(t\right)\right),\cdots ,{\tilde{\alpha }}_R\left({{\mathit{\beta }}_Q}^R\left(t\right)\right)\right]\in {\mathbb{C}}^{N_{RF}^R\times Q}$$

(34)

$${\mathit{U}}_2=\left[{\tilde{\kappa }}_1\left(t\right),\cdots {\tilde{\kappa }}_Q\left(t\right)\right]\in {\mathbb{C}}^{N_{RF}^T\times Q}$$

(35)

$${\mathit{U}}_3=\left[{\tilde{\mathit{\alpha }}}_T\left({{\gamma }_1}^T\left(t\right)\right),\cdots {\tilde{\mathit{\alpha }}}_T\left({{\gamma }_Q}^T\left(t\right)\right)\right]\in {\mathbb{C}}^{\mathrm{T}\times Q}$$

(36)

The FCTPM based on tensor decomposition is utilized the PM (Power Method) for the auxiliary unfolding matrix decomposition. FCTPM is based on a combination of Tensor-train decomposition and power. The FCTPM approximates the most dominant rank-negative-one tensor of the input tensor by the Tensor-train decomposition, and each unfolding matrix is also decomposed using a rank-negative-one matrix decomposition. This property allows each unfolding matrix to ignore the concatenation between other unfolding matrices. Therefore, a cyclic update method is proposed to consider the connectivity between the unfolding matrices. After all the left singular vectors have been computed, the order of the update steps for the right singular vectors needs to be changed. The FCTPM utilizes the PM for the auxiliary unfolding matrix decomposition. Therefore, the PM ensures the convergence of the decomposition factors. The cyclic updating method enhances the connectivity between the decomposition factors.

The FCTPM then consists of two main components, namely the rank-negative-one tensor approximated by the Tensor-train decomposition of the auxiliary expansion matrix, and the PM of determining the left and right singular vectors obtained from the randomly initialized vectors. The core tensor in the Tensor-train decomposition ${\mathcal{G}}^{\left(n\right)}$ calculation is based on the N-dimensional input tensor $\mathcal{Y}\in {\mathbb{R}}^{I_1\times \cdots I_N}$ the low-rank approximation of the auxiliary matrix of $R_N$ denotes the received signal tensor $\mathcal{Y}$ the Tensor-train rank in Tensor-train decomposition. In FCTPM, a rank-negative-one approximation is made to the high-dimensional input tensor using the Tensor-train decomposition. The received signal at theuser is compressed into a low-order core tensor ${\mathcal{G}}^{\left(1\right)}$ , and ${\mathcal{G}}^{\left(2\right)}$ , and ${\mathcal{G}}^{\left(3\right)}$ of the concatenation. The received signal tensor $\mathcal{Y}$ is written in the following form:

$$\mathcal{Y}=\sum _{r_1=1}^{R_1}\cdots \sum _{r_2=1}^{R_2}{\mathcal{G}}_{r_0,:,r_1}^{\left(1\right)}\circ {\mathcal{G}}_{r_1,:,r_2}^{\left(2\right)}\circ {\mathcal{G}}_{r_2,:,r_3}^{\left(3\right)}+\mathcal{N}$$

$$={\mathcal{G}}_{1,:,1}^{\left(1\right)}\circ {\mathcal{G}}_{1,:,1}^{\left(2\right)}\circ {\mathcal{G}}_{1,:,1}^{\left(3\right)}+\mathcal{N}$$

$$={\mathit{g}}^{\left(1\right)}\circ {\mathit{g}}^{\left(2\right)}\circ {\mathit{g}}^{\left(3\right)}+\mathcal{N}$$

(37)

Where ${\mathcal{G}}^{\left(n\right)}$ denotes the core tensor of size $R_{n-1}\times {\mathit{I}}_n\times R_n$ of the core tensor of the rank-negative-one Tensor-train with all elements of the Tensor-train rank being 1, then ${\mathcal{G}}^{\left(n\right)}$ size become $1\times {\mathit{I}}_n\times 1$ the core tensor of the rank-negative-one Tensor-train. Furthermore, if the core tensor ${\mathcal{G}}^{\left(n\right)}$ of dimensionality is compressed, the core tensor can be regarded as factored vector data ${\mathit{g}}^{\left(n\right)}ϵ{\mathbb{R}}^{I_{\mathit{n}}}$ that ${\mathit{g}}^{\left(n\right)}$ is the rank-negative-one Tensor-train core tensor ${\mathcal{G}}_{1,:,1}^{\left(n\right)}$ of the first-order factor vector.

To solve for the factor vector ${\mathit{g}}^{\left(n\right)}$ , equation (32) is expressed as the following equation (38):

$$\underset{{\mathit{g}}^{\left(1\right)},{\mathit{g}}^{\left(2\right)},{\mathit{g}}^{\left(3\right)}}{\min }{‖\mathcal{Y}-d{\mathit{g}}^{\left(1\right)}\circ {\mathit{g}}^{\left(2\right)}\circ {\mathit{g}}^{\left(3\right)}‖}^2$$

$$s_{.} t_{.} d>0,{{\mathit{g}}^{\left(n\right)}}^T{\mathit{g}}^{\left(n\right)}=1\left(n=\mathrm{1,2},3\right)$$

(38)

Where $d$ is the principal singular value of the approximate rank-negative-one tensor. According to the FCTPM, the higher order tensor is decomposed in the Tensor-train decomposition by a low-rank decomposition of the auxiliary expansion matrix of the input tensor.

The solution to the factor vector ${\mathit{g}}^{\left(n\right)}$ can be solved in a uniform manner by expressing:

$$\underset{{\mathit{g}}^{\left(N-1\right)},{\mathit{g}}^{\left(N\right)}}{\min }{‖{\mathit{Y}}_{〈N-1〉}-d^{\left(N-1\right)}{\mathit{g}}^{\left(N-1\right)}{{\mathit{g}}^{\left(N\right)}}^T‖}_F^2$$

$$s_{.} t_{.} d>0,{{\mathit{g}}^{\left(N-1\right)}}^{\mathit{T}}{\mathit{g}}^{\left(N-1\right)}=1, {{\mathit{g}}^{\left(N\right)}}^{\mathit{T}}{\mathit{g}}^{\left(N\right)}=1$$

(39)

Where ${\mathit{Y}}_{〈N-1〉}=reshape\left({\mathit{v}}^{\left(N-2\right)},\left[I_{N-1},I_N\right]\right)，d^{\left(N-1\right)}$ is the principal singular value after reshaping the matrix. The principal singular value of the rank-negative-one tensor can be calculated as $d=\prod _{i=1}^{N-1}{d}^{\left(i\right)}$ that is the corresponding factor vector. ${\mathit{g}}^{\left(N-1\right)}$ and ${\mathit{g}}^{\left(N\right)}$ is the corresponding factor vector, from which the factor vector can then be composed to derive the corresponding factor matrix. The solution process is repeated until the right singular vector in the decomposition process ${\mathit{v}}^{\left(N-1\right)}$ is the last factor vector ${\mathit{g}}^{\left(N\right)}$. Therefore, the first factor vector ${\mathit{g}}^{\left(1\right)}$ can be expressed as follows:

$$\underset{{\mathit{g}}^{\left(1\right)},\mathit{v}}{\min }{‖{\mathit{Y}}_{〈1〉}-d^{\left(1\right)}{\mathit{g}}^{\left(1\right)}{{\mathit{v}}^{\left(1\right)}}^{\mathit{T}}‖}_F^2$$

$$s_{.} t_{.} d^{\left(1\right)}>0,{{\mathit{g}}^{\left(1\right)}}^T{\mathit{g}}^{\left(1\right)}=1, {{\mathit{v}}^{\left(1\right)}}^T{\mathit{v}}^{\left(1\right)}=1$$

(40)

Where ${\mathit{Y}}_{〈1〉}$ is the modulo 1 expansion matrix, $d^{\left(1\right)}$ is the ${\mathit{Y}}_{〈1〉}$ the principal singular values of ${\mathit{g}}^{\left(1\right)}$ represents the first factor vector of the rank-negative-one Tensor-train decomposition, and ${\mathit{v}}^{\left(1\right)}$ is the right singular vector of ${\mathit{Y}}_{〈1〉}$. After obtaining ${\mathit{g}}^{\left(1\right)}$, ${\mathit{v}}^{\left(1\right)}$ is expressed as the modulo 2 expansion matrix form of the received signal tensor $\mathcal{Y}$. The second factor vector ${\mathit{g}}^{\left(2\right)}$ is calculated as:

$${\mathit{Y}}_{〈2〉}=reshape\left({\mathit{v}}^{\left(1\right)},\left[I_2,\prod _{j=3}^N{I}_{\mathit{j}}\right]\right)$$

(41)

${\mathit{Y}}_{〈2〉}$ the rank-negative-one matrix decomposition is expressed as the following equation:

$$\underset{{\mathit{g}}^{\left(2\right)},{\mathit{v}}^{\left(2\right)}}{\min }{‖{\mathit{Y}}_{〈2〉}-d^{\left(2\right)}{\mathit{g}}^{\left(2\right)}{{\mathit{v}}^{\left(2\right)}}^{\mathit{T}}‖}_F^2$$

$$s_{.} t_{.} d^{\left(2\right)}>0,{{\mathit{g}}^{\left(2\right)}}^T{\mathit{g}}^{\left(2\right)}=1, {{\mathit{v}}^{\left(2\right)}}^T{\mathit{v}}^{\left(2\right)}=1$$

(42)

Where $d^{\left(2\right)}$ is the principal singular value of the tensor $\mathcal{Y}$ modulo 2 matrix ${\mathit{Y}}_{〈2〉}$, and ${\mathit{g}}^{\left(2\right)}$ represents the second factor vector of the rank-negative-one Tensor-train decomposition. Similarly, the third factor vector can be obtained. This computation can be performed repeatedly until the right singularity vector of the decomposition process is the last factor vector. To ensure convergence of the factorization, a matrix factorization method PM is incorporated. In solving the first factor matrix, the PM method initializes the factor vector ${\mathit{g}}^{\left(1\right)}$ to a random unit vector satisfying ${‖{\mathit{g}}^{\left(1\right)}‖}_2=1$. After that, the vector matrix operation is performed iteratively until the factor vectors satisfy the pre-set specific termination conditions, as shown below:

$${\mathit{v}}^{\left(1\right)}=\frac{{\mathit{Y}}_{〈1〉}^T{\mathit{g}}^{\left(1\right)}}{{‖{\mathit{Y}}_{〈1〉}^T{\mathit{g}}^{\left(1\right)}‖}_2}ϵ{\mathbb{R}}^{\prod _{j=2}^N{I}_{\mathit{j}}}$$

(43)

$${\mathit{g}}^{\left(1\right)}=\frac{{\mathit{Y}}_{〈1〉}{\mathit{v}}^{\left(1\right)}}{{‖{\mathit{Y}}_{〈1〉}{\mathit{v}}^{\left(1\right)}‖}_2}ϵ{\mathbb{R}}^{I_1}$$

(44)

Where ${\mathit{g}}^{\left(1\right)}$ and ${\mathit{v}}^{\left(1\right)}$ are ${\mathit{Y}}_{〈1〉}$ the left and right singular vectors, respectively. The main singular values $d^{\left(1\right)}$ associated with it can be expressed as $d^{\left(1\right)}={{\mathit{g}}^{\left(1\right)}}^T{\mathit{Y}}_{〈1〉}{\mathit{v}}^{\left(1\right)}$. The left singular vector calculation follows immediately after the right singular vector calculation. Assuming the existence of an orthogonal basis $\overline{\mathit{U}}ϵ{\mathbb{R}}^{I_1\times R}$ and $\overline{\mathit{V}}ϵ{\mathbb{R}}^{I_2\cdots I_N\times R}$, the unfolding matrix ${\mathit{Y}}_{〈1〉}$ satisfies the following equation:

$${\mathit{Y}}_{〈1〉}=\sum _{j=1}^R{d}_j{\overline{\mathit{U}}}_{:,j}{\overline{\mathit{V}}}_{:,j}^T$$

(45)

Where $d_j$ is ${\mathit{Y}}_{〈1〉}$ the singular values in corresponding to the singular vectors. After several iterations, the ${\mathit{v}}^{\left(1\right)}$ and ${\mathit{g}}^{\left(1\right)}$ can be defined as the following equation, respectively:

$${\mathit{v}}^{\left(1\right)}={\delta }_{{\mathit{v}}^{\left(1\right)}}\sum _{j=1}^R{d}_j^{2k+1}{\overline{\mathit{V}}}_{:,j}{\overline{\mathit{U}}}_{:,j}^T{\mathit{g}}^{\left(1\right)}$$

(46)

$${\mathit{g}}^{\left(1\right)}={\delta }_{{\mathit{g}}^{\left(1\right)}}\sum _{j=1}^R{d}_j^{2k}{\overline{\mathit{U}}}_{:,j}{\overline{\mathit{U}}}_{:,j}^T{\mathit{g}}^{\left(1\right)}$$

(47)

Where ${\delta }_{{\mathit{g}}^{\left(1\right)}}$ and${ \delta }_{{\mathit{v}}^{\left(1\right)}}$ are the corresponding normalisation factors. Thus, after several iterations of PM, ${\mathit{v}}^{\left(1\right)}$ and ${\mathit{g}}^{\left(1\right)}$ satisfy $\frac{{‖{\mathit{v}}^{\left(1\right)}‖}_2^2}{{‖{\mathit{g}}^{\left(1\right)}‖}_2^2}=1$ by the convergence property, the $\frac{{\delta }_{{\mathit{v}}^{\left(1\right)}}}{{\delta }_{{\mathit{g}}^{\left(1\right)}}}$ converges to the maximum singular value $d_1$, at ${\mathit{Y}}_{〈1〉}^T{\mathit{g}}^{\left(1\right)}=\left(\frac{{\delta }_{{\mathit{v}}^{\left(1\right)}}}{{\delta }_{{\mathit{g}}^{\left(1\right)}}}\right) {\mathit{v}}^{\left(1\right)}$, ${‖{\mathit{Y}}_{〈1〉}^T{\mathit{g}}^{\left(1\right)}-d_1{\mathit{v}}^{\left(1\right)}‖}_2=0$.

From the previous rank-negative-one Tensor-train decomposition, the nth factor vector ${\mathit{g}}^{\left(n\right)}$ depends on the previously computed right singular vector ${\mathit{v}}^{\left(n-1\right)}$, while ${\mathit{v}}^{\left(n-1\right)}$ is obtained from n-1 factor vector ${\mathit{g}}^{\left(n-1\right)}$computed. Because the matrix decomposition of each factor vector is performed independently. These properties can lead to problems with local optimization, resulting in a loss of critical performance in FCTPM. To enhance the connectivity between the decomposition factors, a circular update method is proposed, which changes the order in which the left singular vectors are computed is sufficient. According to the decomposition method of Tensor-train of rank-negative-one and the nature of PM, ${\mathit{v}}^{\left(n\right)}$ can be obtained from ${\mathit{g}}^{\left(n+1\right)}$ and ${\mathit{v}}^{\left(n+1\right)}$ which are reshaped to form ${\mathit{v}}^{\left(n\right)}$. The relationship ${\mathit{g}}^{\left(n\right)}$ is updated by the following way:

$${\mathit{v}}^{\left(n\right)}\approx {\widehat{\mathit{v}}}^{\left(n\right)}=reshape\left(d^{\left(n+1\right)}{\mathit{g}}^{\left(n+1\right)}{{\mathit{v}}^{\left(n+1\right)}}^T,\left[I_{n+1}\cdots I_N\right]\right)$$

$${\mathit{g}}^{\left(n\right)}=\frac{{\mathcal{Y}}_{〈n〉}{\widehat{\mathit{v}}}^{\left(n\right)}}{{‖{\mathcal{Y}}_{〈n〉}{\widehat{\mathit{v}}}^{\left(n\right)}‖}_2}$$

(48)

The updated factor vector ${\mathit{g}}^{\left(2\right)}$ can be derived as ${\mathit{g}}^{\left(2\right)}=\frac{{\mathcal{Y}}_{〈2〉}{\widehat{\mathit{v}}}^{\left(2\right)}}{{‖{\mathcal{Y}}_{〈2〉}{\widehat{\mathit{v}}}^{\left(2\right)}‖}_2}$ .

Channel parameter is estimated for UAV with millimeter wave massive MIMO based on the FCTPM. The number of factor matrix time slots is associated with the structure of the transmission frame, using a downlink transmission frame structure. It consists of an AOA/AOD estimation phase ($T_1$ time slot), a fast-fading factor estimation phase ($T_2$ time slot) and the data transmission phase ($T_3$ time slot). Using equation (32) and equation (33), the received signal at the terrestrial user is represented as the following equation:

$${\mathcal{Y}}^{T_2}=\mathcal{T};{\mathit{U}}_1,{\mathit{U}}_2,{\mathit{U}}_3+{\mathcal{N}}^{T_2}$$

(49)

Matrix $\left\{{\mathit{U}}_1,{\mathit{U}}_2,{\mathit{U}}_3\right\}$ is the received signal at the user. $\mathcal{Y}$ is corresponding to the matrix of three factors. ${\mathit{U}}_2$ is the corresponding fast fading factor matrix. From the Tensor-train decomposition's$Kruskal$ representation, ${\mathit{U}}_n$ and ${\mathcal{G}}_n$ the correspondence can be expressed by the following equation:

$${\mathit{U}}_n={\left[{\mathcal{G}}_n\right]}_{\left(2\right)}\left(1_{R_{>\left(n\right)}}^T⨂{\mathsf{{\rm I}}}_{R_{n-1}}{R}_n⨂1_{R_{<n-1}}^T\right),n=\mathrm{1,2},3$$

$$={\left[{\mathit{g}}^{\left(n\right)}\right]}_{\left(2\right)}\left(1_{R_{>\left(n\right)}}^T⨂{\mathsf{{\rm I}}}_{R_{n-1}}{R}_n⨂1_{R_{<n-1}}^T\right),n=\mathrm{1,2},3$$

(50)

The factor matrix ${\mathit{U}}_2={\left[{\mathit{g}}^{\left(n\right)}\right]}_{\left(2\right)}\left(1_{R_{>\left(2\right)}}^T⨂{\mathsf{{\rm I}}}_{R_1}{R}_n⨂1_{R_1}^T\right)$. From equation (35), we can find that the factor matrix ${\mathit{U}}_2$ has the fast fading channel power fading factor associated with the Doppler shift in its elements. We can derive the fast fading channel power factor. We can obtain${\mathit{U}}_2=\left[{\tilde{\kappa }}_1,\cdots {\tilde{\kappa }}_Q\right]ϵ{\mathbb{C}}^{T_2\times Q}$, ${\tilde{\kappa }}_Q={\left[{\tilde{\kappa }}_1\left(1\right),\cdots {\tilde{\kappa }}_Q\left(T_2\right)\right]}^Tϵ{\mathbb{C}}^{T_2\times 1}$. According to the correlation scheme, the Doppler shift on the $q$ path is solved as follows:

$${\widehat{f}}_q=\underset{f_{\mathrm{q}}}{\mathrm{argmax}}\frac{\left|[{\widehat{U}}_2{]}_{:,q}^H\kappa \left(f_{\mathrm{q}}\right)\right|}{{‖[{\widehat{U}}_2{]}_{:,q}‖}_2{‖\kappa \left(f_{\mathrm{q}}\right)‖}_2}.1\le \mathrm{q}\le Q . {\widehat{f}}_q\in [0,f_{max}]$$

(51)

$$f_{max}=\frac{f_c{v}^{\mathrm{\text{'}}}}{c}{\kappa }_q\left(t\right)=b_q{e}^{j2\pi f_{q}t}$$

(52)

Where $f_{max}$ is the maximum Doppler shift,${ f}_c$ is the carrier frequency, $v^{\mathrm{\text{'}}}$ is the speed of the drone, and $c$ is the speed of light. The path gain is then estimated from the estimated Doppler shift, i.e:

$${\widehat{b}}_q=\left[\kappa \left({\widehat{f}}_q\right){]}^{\uparrow }\right[{\widehat{U}}_2{]}_{:,q}$$

(53)

Finally, all parameters of the time-varying channel of the millimeter-wave MIMO in UAV are estimated and the channel matrix can be recovered from (51).

4. Hybrid beamforming and RIS phase shift matrix design

4.1. Hybrid beamforming and RIS phase shift matrix design based on Tensor-train decomposition

Define$\sqrt{\frac{{{\mathrm{T}}^{\mathrm{‘}}N}_R}{Q_{{\ddot{M}}_k}}}$ and$\sqrt{\frac{{{\mathrm{T}}^{\mathrm{‘}}N}_T}{Q_{\ddot{R}}}}$ as constants and set it to 1. The millimeter-wave massive channel from the RIS to the user side ${\ddot{\mathit{M}}}_k$ and the millimeter wave massive channel from the UAV side to the RIS $\ddot{\mathit{R}}$ can be written with Tensor-train decomposition form as follows, respectively :

$${\ddot{\mathcal{M}}}_k=\sum _{q=1}^{Q_{{\ddot{M}}_k}}\mathcal{T}{\stackrel{̿}{\kappa }}_{km,q}·{\stackrel{̿}{\mathit{\alpha }}}_R\left({{\phi }_{km,q}}^R\left(t\right)\right)·{\stackrel{̿}{\mathit{\alpha }}}_{P_1}^H\left({{\phi }_{km,q}}^T\left(t\right),{{\eta }_{km,q}}^T\left(t\right)\right)$$

$$=\mathcal{T};{\stackrel{̿}{\mathit{A}}}_{R,km},{\stackrel{̿}{\mathsf{\Lambda }}}_{km,q},{\stackrel{̿}{\mathit{A}}}_{P_1,km}^H$$

$$={\stackrel{̿}{\mathit{g}}}_{R,km}^{\left(1\right)}\circ {\stackrel{̿}{\mathit{g}}}_{km,q}^{\left(2\right)}\circ {\stackrel{̿}{\mathit{g}}}_{P_1,km}^{H\left(3\right)}$$

(54)

$$\ddot{\mathcal{R}}=\sum _{q=1}^{Q_{\ddot{R}}}{\mathcal{T}\stackrel{̿}{\kappa }}_{r,q}\left(t\right)·{\stackrel{̿}{\mathit{\alpha }}}_{P_2}\left({{\phi }_{r,q}}^R\left(t\right),{{\eta }_{r,q}}^R\left(t\right)\right)·{\stackrel{̿}{\mathit{\alpha }}}_T^H\left({{\nu }_{r,q}}^T\left(t\right)\right)$$

$$=\mathcal{T};{\stackrel{̿}{\mathit{A}}}_{P_2,r},{\stackrel{̿}{\mathsf{\Lambda }}}_{r,q},{\stackrel{̿}{\mathit{A}}}_{T,r}^H$$

$$={\stackrel{̿}{\mathit{g}}}_{P_2,r}^{\left(1\right)}\circ {\stackrel{̿}{\mathit{g}}}_{r,q}^{\left(2\right)}\circ {\stackrel{̿}{\mathit{g}}}_{T,r}^{H\left(3\right)}$$

(55)

Where ${\stackrel{̿}{\kappa }}_{km,q}\left(t\right)={\left[{\kappa }_{km,q}\left(1\right),\cdots ,{\kappa }_{km,q}\left(T\right)\right]}^Tϵ{\mathbb{C}}^{T\times 1}$, ${\stackrel{̿}{\kappa }}_{r,q}\left(t\right)={\left[{\kappa }_{r,q}\left(1\right),\cdots ,{\kappa }_{r,q}\left(T\right)\right]}^Tϵ{\mathbb{C}}^{T\times 1}$. The three factor matrices obtained by tensor decomposition of ${\ddot{\mathit{M}}}_k$ are the receive antenna array response matrix ${\stackrel{̿}{\mathit{A}}}_{R,km}$, the path gain factor matrix ${\mathsf{\Lambda }}_{km,q}$, the RIS reflector antenna array response matrix ${\stackrel{̿}{\mathit{A}}}_{P_1,km}^H$. The three factor matrices are defined as follows, respectively :

$${\stackrel{̿}{\mathit{A}}}_{R,km}={\stackrel{̿}{\mathit{g}}}_{R,km}^{\left(1\right)}\left(1_{R_{>\left(2\right)}}^T⨂{\mathsf{{\rm I}}}_{R_1}{R}_{R,km}⨂1_{R_1}^T\right)$$

$$=\left[{\stackrel{̿}{\mathit{\alpha }}}_R（{{\phi }_{km,1}}^R\left(t\right)）,\cdots ,{\stackrel{̿}{\mathit{\alpha }}}_R（{{\phi }_{km,Q_{{\ddot{M}}_k}}}^R\left(t\right)）\right]ϵ{\mathbb{C}}^{N_R\times Q_{{\ddot{M}}_k}}$$

(56)

$${\mathsf{\Lambda }}_{km,q}={\stackrel{̿}{\mathit{g}}}_{km,q}^{\left(2\right)}\left(1_{R_{>\left(2\right)}}^T⨂{\mathsf{{\rm I}}}_{R_1}{R}_{km,q}⨂1_{R_1}^T\right)$$

$$=\left[{\stackrel{̿}{\kappa }}_{km,1}\left(t\right),\cdots ,{\stackrel{̿}{\kappa }}_{km,Q_{{\ddot{M}}_k}}\left(t\right)\right]ϵ{\mathbb{C}}^{T\times Q_{{\ddot{M}}_k}}$$

(57)

$${\stackrel{̿}{\mathit{A}}}_{P_1,km}^H={\stackrel{̿}{\mathit{g}}}_{P_1,km}^{H\left(3\right)}\left(1_{R_{>\left(2\right)}}^T⨂{\mathsf{{\rm I}}}_{R_1}{R}_{km,q}⨂1_{R_1}^T\right)$$

$$=\left[{\stackrel{̿}{\mathit{\alpha }}}_{P_1}\left({{\phi }_{km,1}}^T\left(t\right),{{\eta }_{km,1}}^T\left(t\right)\right),\cdots ,{\stackrel{̿}{\mathit{\alpha }}}_{P_1}\left({{\phi }_{km,Q_{{\ddot{M}}_k}}}^T\left(t\right),{{\eta }_{km,Q_{{\ddot{M}}_k}}}^T\left(t\right)\right)\right]ϵ{\mathbb{C}}^{{\mathrm{T}}^{\mathrm{‘}}\times Q_{{\ddot{M}}_k}}$$

(58)

Where ${\ddot{\mathit{R}}}_k$ is also obtained by tensor decomposition to obtain the RIS receive antenna array response factor matrix ${\stackrel{̿}{\mathit{A}}}_{P,r}$, the fast-fading factor matrix ${\stackrel{̿}{\mathsf{\Lambda }}}_{r,q}$, and the UAV end transmit antenna factor matrix ${\stackrel{̿}{\mathit{A}}}_{T,r}^H$. The three factor matrices are defined respectively as follows:

$${\stackrel{̿}{\mathit{A}}}_{P,r}={\stackrel{̿}{\mathit{g}}}_{P,r}^{\left(1\right)}\left(1_{R_{>\left(2\right)}}^T⨂{\mathsf{{\rm I}}}_{R_1}{R}_{P,r}⨂1_{R_1}^T\right)$$

$$=\left[{\stackrel{̿}{\mathit{\alpha }}}_P\left({{\phi }_{r,1}}^R\left(t\right),{{\eta }_{r,1}}^R\left(t\right)\right),\cdots ,{\stackrel{̿}{\mathit{\alpha }}}_P\left({{\phi }_{r,Q_{\ddot{R}}}}^R\left(t\right),{{\eta }_{r,Q_{\ddot{R}}}}^R\left(t\right)\right)\right]ϵ{\mathbb{C}}^{{\mathrm{T}}^{\mathrm{‘}}\times Q_{\ddot{R}}}$$

(59)

$${\stackrel{̿}{\mathsf{\Lambda }}}_{r,q}={\stackrel{̿}{\mathit{g}}}_{r,q}^{\left(2\right)}\left(1_{R_{>\left(2\right)}}^T⨂{\mathsf{{\rm I}}}_{R_1}{R}_{r,q}⨂1_{R_1}^T\right)$$

$$=\left[{\stackrel{̿}{\kappa }}_{km,1}\left(t\right),\cdots ,{\stackrel{̿}{\kappa }}_{km,Q_{\ddot{R}}}\left(t\right)\right]ϵ{\mathbb{C}}^{T\times Q_{\ddot{R}}}$$

(60)

$${\stackrel{̿}{\mathit{A}}}_{T,r}^H={\stackrel{̿}{\mathit{g}}}_{T,r}^{H\left(3\right)}\left(1_{R_{>\left(2\right)}}^T⨂{\mathsf{{\rm I}}}_{R_1}{R}_{T,r}⨂1_{R_1}^T\right)$$

$$=\left[{\stackrel{̿}{\mathit{\alpha }}}_P\left({{\nu }_{r,1}}^T\left(t\right)\right),\cdots ,{\stackrel{̿}{\mathit{\alpha }}}_P\left({{\nu }_{r,Q_{\ddot{R}}}}^T\left(t\right)\right)\right] ϵ{\mathbb{C}}^{N_T\times Q_{\ddot{R}}}$$

(61)

According to the existing tensor form ${\ddot{\mathcal{M}}}_k$ and tensor form $\ddot{\mathcal{R}}$, the received signal ${\mathcal{Y}}_k$ The tensor form of k is expressed as:

$${\mathcal{Y}}_k={\mathit{W}}_{BB,k}^H{\mathit{W}}_{RF,k}^H·{\ddot{\mathcal{M}}}_k\mathbf{\Phi }\ddot{\mathcal{R}}·{\mathit{F}}_{RF}^T{\mathit{F}}_{BB,k}^T{\mathit{s}}_k+\sum _{j\ne k}^K{\mathit{W}}_{BB,k}^H{\mathit{W}}_{RF,k}^H·{\ddot{\mathcal{M}}}_k\mathbf{\Phi }\ddot{\mathcal{R}}{·\mathit{W}}_{BB,k}^H{\mathit{W}}_{RF,k}^H·{\mathit{F}}_{RF}^T{\mathit{F}}_{BB,j}^T{\mathit{s}}_j+{\mathcal{N}}_k$$

$$=\mathcal{T};{\mathit{W}}_{BB,k}^H{\mathit{W}}_{RF,k}^H,{\stackrel{̿}{\mathit{A}}}_{R,km},{\stackrel{̿}{\mathsf{\Lambda }}}_{km,q},{\stackrel{̿}{\mathit{A}}}_{P,km}^H,\mathbf{\Phi },\mathcal{T};{\stackrel{̿}{\mathit{A}}}_{P,r},{\stackrel{̿}{\mathsf{\Lambda }}}_{r,q},{\stackrel{̿}{\mathit{A}}}_{T,r}^H,{\mathit{F}}_{RF}^T{\mathit{F}}_{BB,k}^T{\mathit{s}}_k+$$

$${\displaystyle {\sum }_{j\ne k}^{\mathit{K}}}\mathcal{T};{\mathit{W}}_{BB,k}^H{\mathit{W}}_{RF,k}^H,{\stackrel{̿}{\mathit{A}}}_{R,km},{\stackrel{̿}{\mathsf{\Lambda }}}_{km,q},{\stackrel{̿}{\mathit{A}}}_{P,km}^H,\mathbf{\Phi },\mathcal{T};{\stackrel{̿}{\mathit{A}}}_{P,r},{\stackrel{̿}{\mathsf{\Lambda }}}_{r,q},{\stackrel{̿}{\mathit{A}}}_{T,r}^H,{\mathit{F}}_{RF}^T{\mathit{F}}_{BB,j}^T{\mathit{s}}_j+{\mathcal{N}}_k$$

(62)

Where ${\mathcal{Y}}_k\in {\mathbb{C}}^{N_s\times k{N}_{s,k}\times T}$ is the received signal tensor of the kth terrestrial user, and ${\mathcal{N}}_k\in {\mathbb{C}}^{N_s\times k{N}_{s,k}\times T}$ is the noise tensor of the kth terrestrial user.

By utilizing RIS assisted millimeter wave massive MIMO of UAV communication, the spectral efficiency $R_k$ of hybrid beamforming can be expressed as:

$$R_k={\displaystyle {\sum }_{k=1}^K}{\log }_2\left|{\mathit{I}}_{k{N}_s}+\frac{p^{\mathrm{‘}}}{k{N}_s}{\mathit{R}}_{kn}^{-1}{\mathit{W}}_{BB,k}^H{\mathit{W}}_{RF,k}^H{\ddot{\mathit{M}}}_k\mathbf{\Phi }{{\ddot{\mathit{R}}}_k\mathit{F}}_{RF}{\mathit{F}}_{BB,k}{\left({\mathit{W}}_{BB,k}^H{\mathit{W}}_{RF,k}^H{\ddot{\mathit{M}}}_k\mathbf{\Phi }{\ddot{\mathit{R}}}_k{\mathit{F}}_{RF}{\mathit{F}}_{BB,k}\right)}^H\right|$$

(63)

Where ${\mathit{R}}_{kn}={\sigma }_n^2{\mathit{W}}_{BB,k}^H{\mathit{W}}_{RF,k}^H{\mathit{W}}_{RF,k}{\mathit{W}}_{BB,k}$ denotes the $k$ noise covariance matrix at the received user, and $p^{\mathrm{‘}}$ denotes the transmit power. So the maximum spectral efficiency problem can be expressed as the following equation:

$$\underset{{\mathit{F}}_{\mathit{RF}},{\mathit{F}}_{\mathit{BB},k},{\mathit{W}}_{\mathit{RF},k}{,\mathit{W}}_{\mathit{BB},k},\mathbf{\Phi }}{\max }{R}_k$$

$$s_{.}{t_{.}‖{\mathit{F}}_{RF}{\mathit{F}}_{BB,k}‖}_F^2={kN}_{s,k}$$

$${\theta }_{t^{\mathrm{‘}}}ϵ\mathcal{I},\forall t^{\mathrm{‘}}$$

(64)

This is a non-convex NP-hard problem. It decomposes the original optimization problem into two independent optimization subproblems, namely $\mathbf{\Phi }$ the optimization problem and (${\mathit{F}}_{RF},{\mathit{F}}_{BB,k},{\mathit{W}}_{RF,k}{\mathit{W}}_{BB,k}）$ of the optimization problem. Setting up the hybrid channel ${\mathit{H}}_{kc}={\ddot{\mathit{M}}}_k\mathbf{\Phi }{\ddot{\mathit{R}}}_k$. Next, it is necessary to calibrate the RIS phase shift matrix $\mathbf{\Phi }$ Perform optimization solutions.

4.2. RIS phase shift matrix design

To simplify the problem, we consider the use of a fully digital beamforming matrix ${\mathit{F}}_k^fϵ{\mathbb{C}}^{N_T\times N_s}$ and a fully digital combination matrix ${\mathit{W}}_k^fϵ{\mathbb{C}}^{N_R\times N_s}$ to replace the hybrid beamforming matrix ${\mathit{F}}_{RF}$ , the ${\mathit{F}}_{BB,k}$ and the hybrid combiner matrix ${\mathit{W}}_{RF,k}{\mathit{W}}_{BB,k}$ . The previous optimization problem can be expressed as follows:

$$\underset{\mathit{F},\mathit{W},\mathbf{\Phi }}{\max }{\log }_{2}det\left({\mathit{I}}_{k{N}_s}+\frac{{\rho }^{\mathrm{‘}}}{k{N}_{\mathit{s}}{\sigma }_{\mathit{n}}^2}{{\mathit{W}}_k^f}^{-1}{{\ddot{\mathit{M}}}_k}^H\mathbf{\Phi }\times {\ddot{\mathit{R}}}_k{\mathit{F}}_k^f{{\mathit{F}}_k^f}^H{{\ddot{\mathit{R}}}_k}^H{\mathbf{\Phi }}^H{\ddot{\mathit{M}}}_k{\mathit{W}}_k^f\right)$$

$$s_{.}{t_{.}‖{\mathit{F}}_k^f‖}_F^2={kN}_{s,k}$$

$${\theta }_{t^{\mathrm{‘}}}ϵ\mathcal{I}, \forall t^{\mathrm{‘}}$$

(65)

RIS reflection matrix $\mathbf{\Phi }$ is set as fixed. The mixed channel matrix is constructed as ${\mathit{H}}_{kc}={\ddot{\mathit{M}}}_k\mathbf{\Phi }{\ddot{\mathit{R}}}_k$. By decomposing the effective mixed channel matrix ${\mathit{H}}_{kc}$ performing a singular value decomposition, the optimal ${\mathit{F}}_k^f$ and ${\mathit{W}}_k^f$ can be obtained. ${\mathit{H}}_{kc}={\overbrace{\mathit{Q}}}_{kc}{\tilde{\mathsf{\Sigma }}}_{kc}{{\stackrel{̿}{\mathsf{\Gamma }}}_{kc}^{\mathrm{\text{'}}}}^H$. Where ${\overbrace{\mathit{Q}}}_{kc}ϵ{\mathbb{C}}^{N_R\times N_R}$ , ${\stackrel{̿}{\mathsf{\Gamma }}}_{kc}^{\mathrm{\text{'}}}ϵ{\mathbb{C}}^{N_T\times N_T}$ , and ${\tilde{\mathsf{\Sigma }}}_{kc}$ is a rectangular diagonal matrix composed of singular values in descending order. Due to the sparsity of the UAV channel for millimeter wave massive MIMO, the mixed channel matrix ${\mathit{H}}_{kc}$ usually has a low rank. ${\mathit{H}}_{kc}$ approximation can be obtained as ${\mathit{H}}_{kc}\approx {\overbrace{\mathit{Q}}}_{kc1}{\tilde{\mathsf{\Sigma }}}_{kc1}{{\stackrel{̿}{\mathsf{\Gamma }}}_{kc1}^{\text{'}}}^H$. It simultaneously satisfies ${\overbrace{\mathit{Q}}}_{kc1}\triangleq {\overbrace{\mathit{Q}}}_{kc}\left(:,1:N_{s,k}\right)$, ${\tilde{\mathsf{\Sigma }}}_{kc1}\triangleq \tilde{\mathsf{\Sigma }}\left(1:N_{s,k},1:N_{s,k}\right)$ and ${\stackrel{̿}{\mathsf{\Gamma }}}_{kc1}^{\text{'}}\triangleq {\stackrel{̿}{\mathsf{\Gamma }}}_{kc}^{\text{'}}\left(:,1:N_{s,k}\right)$. Afterwards, with a fixed RIS reflection matrix, the channel matrix ${\mathit{H}}_{kc}$ of the optimal unconstrained beamforming matrix is ${{\mathit{F}}_k^f}_{opt}={\stackrel{̿}{\mathsf{\Gamma }}}_{kc1}^{\text{'}}$, ${{\mathit{W}}_k^f}_{opt}$ =${\overbrace{\mathit{Q}}}_{kc1}$. Substituting this into equation (65), a simplified optimization problem can be obtained as follows:

$$\underset{\mathbf{\Phi }}{\max }{\log }_{2}det\left({\mathit{I}}_{k{N}_s}+\frac{{\rho }^{\mathrm{‘}}}{k{N}_{s,k}{\sigma }_n^2}{\tilde{\mathsf{\Sigma }}}_{\mathrm{k}\mathrm{c}1}^2\right)$$

$$s_{.}{t}_{.}{\theta }_{t^{\mathrm{‘}}}ϵ\mathcal{I}, \forall t^{\mathrm{‘}}$$

(66)

This optimization problem remains intractable due to the determinant operation of the objective function and the non-convex constraint of the RIS reflection. To simplify the problem, an upper bound on the above objective is first derived to eliminate the determinant function as follows:

$${\log }_{2}det\left({\mathit{I}}_{k{N}_s}+\frac{{\rho }^{\mathrm{‘}}}{k{N}_{s,k}{\sigma }_n^2}{\tilde{\mathsf{\Sigma }}}_{\mathrm{k}\mathrm{c}1}^2\right)\begin{array}{c}\left(a\right)\\ \le \end{array}{N}_{s,k}{\log }_2\left(1+\frac{{\rho }^{\mathrm{‘}}}{k{N_{s,k}}^2{\sigma }_n^2}\mathrm{T}\mathrm{r}\left({\tilde{\mathsf{\Sigma }}}_{\mathrm{k}\mathrm{c}1}^2\right)\right)$$

$$\begin{array}{c}\left(b\right)\\ \le \end{array}{N}_{s,k}{\log }_2\left(1+\frac{{\rho }^{\mathrm{‘}}}{k{N_{s,k}}^2{\sigma }_n^2}\mathrm{T}\mathrm{r}\left({\mathit{H}}_{kc}{\mathit{H}}_{kc}^H\right)\right)$$

(67)

Where $\left(a\right)$ according to Jensen's inequality holds. When$\left(b\right)$ holds, $\mathrm{T}\mathrm{r}\left({\mathit{H}}_{kc}{\mathit{H}}_{kc}^H\right)={\sum }_{i=1}^{\tilde{Q}}{\lambda }_i^2\ge {\sum }_{i=1}^{N_s}{\lambda }_i^2=\mathrm{T}\mathrm{r}\left({\tilde{\mathsf{\Sigma }}}_{\mathrm{k}\mathrm{c}1}^2\right)$, at t$rank\left({\tilde{\mathsf{\Sigma }}}_{kc}\right)=N_s$ the equal sign holds, ${\tilde{\lambda }}_{kci}$ denotes $\tilde{\mathsf{\Sigma }}$ the first $i$ diagonal element. The reflection design is then used to maximize this upper bound, as in the following equation:

$$\underset{\mathbf{\Phi }}{\max }Tr\left({\mathit{H}}_{kc}{\mathit{H}}_{kc}^H\right)$$

$$=\underset{\mathbf{\Phi }}{\max }Tr\left({\ddot{\mathit{M}}}_k^H\mathbf{\Phi }{\ddot{\mathit{R}}}_k{\ddot{\mathit{R}}}_k^H{\mathbf{\Phi }}^H{\ddot{\mathit{M}}}_k\right)$$

$$\begin{array}{c}\left(a\right)\\ =\end{array}\underset{\mathbf{\Phi }}{\mathrm{m}\mathrm{a}\mathrm{x}}\sum _{i=1}^{N_R}{{{\ddot{\mathit{m}}}_k}_i}^H\mathbf{\Phi }{\ddot{\mathit{R}}}_k{{\ddot{\mathit{R}}}_k}^H{\mathbf{\Phi }}^H{\ddot{\mathit{m}}}_{ki}$$

$$\begin{array}{c}\left(b\right)\\ =\end{array}\underset{\mathbf{\Phi }}{\mathrm{m}\mathrm{a}\mathrm{x}}{\sum }_{i=1}^{N_R}{{\stackrel{̿}{\mathsf{\Gamma }}}_{kc}^{\text{'}}}^H\ddot{\mathit{T}}{\stackrel{̿}{\mathsf{\Gamma }}}_{kc}^{\text{'}}$$

(68)

When $\left(a\right)$ holds, ${{\ddot{\mathit{m}}}_k}_i（\forall i=1,\cdots ,N_R$ ) is the channel matrix ${\ddot{\mathit{M}}}_k$ of the $i$ column. When$\left(b\right)$ holds, the matrix $\ddot{\mathit{T}}$ is defined $\ddot{\mathit{T}}\triangleq \sum _{i=1}^{N_R}diag\left({{{\ddot{\mathit{m}}}_k}_i}^H\right){\ddot{\mathit{R}}}_k{{\ddot{\mathit{R}}}_k}^{H}diag\left({{\ddot{\mathit{m}}}_k}_i\right)$.

To solve this problem efficiently, we use other fixed other phase shifts to optimize each reflection element in turn in an iterative manner until convergence is reached. Specifically, for the $t^{\mathrm{‘}}$ reflection matrix $\mathbf{\Phi }$ element, the corresponding objective function is designed to be linear as follows:

$$2\mathcal{R}\left\{e^{j{\varphi }_{t^{\mathrm{‘}}}}{\mathfrak{I}}_{t^{\mathrm{‘}}}\right\}+{\sum }_{n\ne t^{\mathrm{‘}}}^{{\mathcal{T}}^{\mathrm{‘}}}{\sum }_{i\ne t^{\mathrm{‘}}}^{{\mathcal{T}}^{\mathrm{‘}}}\ddot{\mathit{T}}\left(n,i\right)e^{\mathit{j}\left({\varphi }_n^{\mathrm{\text{'}}}-{\varphi }_i^{\mathrm{\text{'}}}\right)}+\ddot{\mathit{T}}\left(t^{\mathrm{‘}},t^{\mathrm{‘}}\right)$$

(69)

Where ${\mathfrak{I}}_{t^{\mathrm{‘}}}={\sum }_{n\ne t^{\mathrm{‘}}}^{{\mathcal{T}}^{\mathrm{‘}}}\ddot{\mathit{T}}\left(t^{\mathrm{‘}},n\right)e^{-\mathit{j}{\varphi }_{\mathrm{n}}}=\left|{\mathfrak{I}}_{t^{\mathrm{‘}}}\right|e^{-\mathit{j}{\phi }_{t^{\mathrm{‘}}}}$ , the last two terms in the above equation are constants with fixed other reflective elements. The question of the $t^{\mathrm{‘}}$ element can be written as follows:

$$\begin{gathered} \underset{{\varphi }_{t^{\mathrm{‘}}}}{\mathrm{m}\mathrm{a}\mathrm{x}}\mathcal{R}\left\{e^{j{\varphi }_{t^{\mathrm{‘}}}}{\mathfrak{I}}_{t^{\mathrm{‘}}}\right\}=\underset{{\theta }_{t^{\mathrm{‘}}}}{\mathrm{m}\mathrm{a}\mathrm{x}}\mathcal{R}\left\{{\left|{\mathfrak{I}}_{t^{\mathrm{‘}}}\right|e}^{j\left({\varphi }_{t^{\mathrm{‘}}}-{\phi }_{t^{\mathrm{‘}}}\right)}\right\} \\ =\underset{{\theta }_{t^{\mathrm{‘}}}}{\mathrm{m}\mathrm{i}\mathrm{n}}\left|{\varphi }_{t^{\mathrm{‘}}}-{\phi }_{t^{\mathrm{‘}}}\right| \end{gathered}$$

(70)

Where the closed-form solution of the above equation is ${{\varphi }_{t^{\mathrm{‘}}}}^{\mathrm{*}}=\left[\frac{{\phi }_{t^{\mathrm{‘}}}}{∆\theta }\right]\times ∆\varphi$ , $∆\varphi =\raisebox{1ex}{$2\pi $}\!\left/ \!\raisebox{-1ex}{$2^{\mathcal{I}}$}\right.$ is the angular resolution. The phase shift of each reflection element can be designed in turn.

5. Low-complexity hybrid beamforming

5.1. Hybrid beamforming based on spectral efficiency maximization

When the computation of the RIS reflection matrix $\mathbf{\Phi }$ is complete, we will next solve the optimization problem of (${\mathit{F}}_{RF},{\mathit{F}}_{BB,k},{\mathit{W}}_{RF,k}{\mathit{W}}_{BB,k})$. The computational complexity becomes high if the joint optimization is performed at both (${\mathit{F}}_{RF},{\mathit{F}}_{BB,k},{\mathit{W}}_{RF,k}{\mathit{W}}_{BB,k}）$ transceivers at the same time. Because of the constant mode constraints of ${\mathit{F}}_{RF,k}$ and ${\mathit{W}}_{RF,k}$, this optimization problem will be nonconvex, further increasing the complexity of the problem. To facilitate the solution, we decouple the transceiver and the transmitter, while considering the design of the hybrid beamforming matrix at the transmitter side only. The problem in this section is translated into designing the optimal hybrid beamforming matrix $（{\mathit{F}}_{RF},{\mathit{F}}_{BB,k}）.$ This maximizes the system spectral efficiency. After designing the RIS reflection matrix $\mathbf{\Phi }$, the channel ${\mathit{H}}_{kc}={\ddot{\mathit{M}}}_k\mathbf{\Phi }{\ddot{\mathit{R}}}_k$. The channel ${\mathit{H}}_{kc}$ perform a singular value decomposition as the following equation:

$${\mathit{H}}_{kc}={\mathit{G}\text{'}}_{kc}{\mathsf{\Sigma }\mathrm{\text{'}}}_{kc}{\mathit{V}\text{'}}_{kc}^H$$

(71)

Where ${\mathsf{\Sigma }\mathrm{\text{'}}}_{kc}=\left[\begin{array}{cc}{\mathsf{\Sigma }\mathrm{\text{'}}}_{kc,1}& 0\\ 0& {\mathsf{\Sigma }\mathrm{\text{'}}}_{kc,2}\end{array}\right]$ ,${\mathit{V}\text{'}}_{kc}=\left[\begin{array}{cc}{\mathit{V}\text{'}}_{kc,1},& {\mathit{V}\text{'}}_{kc,2}\end{array}\right]$ ,${\mathsf{\Sigma }\mathrm{\text{'}}}_{kc,1}$ represents a stream of information to be processed, and ${\mathit{V}\text{'}}_{kc,1}ϵ{\mathbb{C}}^{N_T\times N_{s,k}}$ represents the product between the information stream to be processed and the number of transmitting antennas. It is assumed that the beamforming matrix received at the receiving end is perfect, that is${\mathit{G}\mathit{\text{'}}}_{kc}={\mathit{W}}_{RF,k}{\mathit{W}}_{BB,k}$ . By referring to the noise covariance matrix of the receiving user above, we can express ${\mathit{R}}_n$ as:

$${\mathit{R}}_n={\sigma }_n^2{\mathit{W}}_{BB,k}^H{\mathit{W}}_{RF,k}^H{\mathit{W}}_{RF,k}{\mathit{W}}_{BB,k}$$

$$={\sigma }_n^2{{\mathit{G}\mathit{\text{'}}}_{kc}}^H{\mathit{G}\mathit{\text{'}}}_{kc}$$

$$={\sigma }_n^2{\mathit{I}}_{{kN}_{s,k}}$$

(72)

Combined with the hybrid channel ${\mathit{H}}_{kc}$, we can express the spectral efficiency of the system through equation (63) as:

$$\begin{gathered} R={\mathit{log}}_2\left(⌈\begin{array}{c}{\mathit{I}}_{k{\mathit{N}}_{s,k}}+\frac{{\rho }^{\mathrm{‘}}}{k{\mathit{N}}_{s,k}}{\mathit{R}}_n^{-1}{\mathit{W}}_{BB,k}^H{\mathit{W}}_{RF,k}^H\left({\mathit{G}\mathbf{\text{'}}}_{kR}{\mathit{\Sigma }\mathbf{\text{'}}}_{kR}{\mathit{V}\mathbf{\text{'}}}_{kc}^H\right)\\ \times {\mathit{F}}_{RF}{\mathit{F}}_{BB,k}{\left({\mathit{W}}_{BB,k}^H{\mathit{W}}_{RF,k}^H\left({\mathit{G}\mathbf{\text{'}}}_{kR}{\mathit{\Sigma }\mathbf{\text{'}}}_{kR}{\mathit{V}\mathbf{\text{'}}}_{kR}^H\right){\mathit{F}}_{RF}{\mathit{F}}_{BB,k}\right)}^H\end{array}⌉\right) \\ ={\mathit{log}}_2\left(⌈{\mathit{I}}_{k{N}_{s,k}}+\frac{{\mathit{\rho }}^{\mathrm{‘}}}{\mathit{k}{\mathit{N}}_{\mathit{s},\mathit{k}}{\mathit{\sigma }}_{\mathit{n}}^2}\left({{\mathit{\Sigma }\mathbf{\text{'}}}_{kR}}^2\right){{\mathit{V}\mathbf{\text{'}}}_{kR}^H\mathit{F}}_{RF}{\mathit{F}}_{BB,k}{{\mathit{F}}_{RF}}^H{{\mathit{F}}_{BB,k}}^H{\mathit{V}\mathbf{\text{'}}}_{kR}⌉\right) \end{gathered}$$

(73)

From the above, it is clear that ${\mathit{V}\mathit{\text{'}}}_{kR,1}ϵ{\mathbb{C}}^{N_T\times N_s}$ . Then the best unconstrained digital beamforming matrix in the system can be ${\mathit{F}}_{opt}={\mathit{V}\mathit{\text{'}}}_{kR,1}$. Because of the constant mode constraint, the ${\mathit{F}}_{opt}$ is not possible to use ${\mathit{F}}_{RF,k}{\mathit{F}}_{BB,k}$ be represented. Assume that ${\mathit{V}\mathit{\text{'}}}_{1c}^H{\mathit{F}}_{RF,k}{\mathit{F}}_{BB,k}\approx {\mathit{I}}_{k{N}_s}$, ${\mathit{V}\mathit{\text{'}}}_{2c}^H{\mathit{F}}_{RF,k}{\mathit{F}}_{BB,k}\approx 0$, the precoding matrix is expressed as:

$${{\mathit{V}\mathbf{\text{'}}}_{kc}^H\mathit{F}}_{RF}{\mathit{F}}_{BB,k}{{\mathit{F}}_{RF}}^H{{\mathit{F}}_{BB,k}}^H{\mathit{V}\mathbf{\text{'}}}_{kc}$$

$$=\left[\begin{array}{cc}{\mathit{V}\mathbf{\text{'}}}_{kc,1}^H{\mathit{F}}_{RF}{\mathit{F}}_{BB,k}{{\mathit{F}}_{RF}}^H{{\mathit{F}}_{BB,k}}^H{\mathit{V}\mathbf{\text{'}}}_{kc,1}& {\mathit{V}\mathbf{\text{'}}}_{kc,1}^H{\mathit{F}}_{RF}{\mathit{F}}_{BB,k}{{\mathit{F}}_{RF}}^H{{\mathit{F}}_{BB,k}}^H{\mathit{V}\mathbf{\text{'}}}_{kc,2}\\ {\mathit{V}\mathbf{\text{'}}}_{kc,2}^H{\mathit{F}}_{RF}{\mathit{F}}_{BB,k}{{\mathit{F}}_{RF}}^H{{\mathit{F}}_{BB,k}}^H{\mathit{V}\mathbf{\text{'}}}_{kc,1}& {\mathit{V}\mathbf{\text{'}}}_{kc,2}^H{\mathit{F}}_{RF}{\mathit{F}}_{BB,k}{{\mathit{F}}_{RF}}^H{{\mathit{F}}_{BB,k}}^H{\mathit{V}\mathbf{\text{'}}}_{kc,2}\end{array}\right]$$

$$\approx \left[\begin{array}{cc}\overline{Q}& 0\\ 0& 0\end{array}\right]$$

(74)

The spectral efficiency can then be approximated by deducing as follows:

$$\begin{gathered} R={\log }_2\left(\left|{\mathit{I}}_{k{N}_{s,k}}+\frac{{\rho }^{\mathrm{‘}}}{k{N}_{s,k}{\sigma }_n^2}\left({{\mathsf{\Sigma }\mathbf{\text{'}}}_{kc}}^2\right){{\mathit{V}\mathbf{\text{'}}}_{kc}^H\mathit{F}}_{RF}{\mathit{F}}_{BB,k}{{\mathit{F}}_{RF}}^H{{\mathit{F}}_{BB,k}}^H{\mathit{V}\mathbf{\text{'}}}_{kc}\right|\right) \\ \approx {\log }_2\left(\left|{\mathit{I}}_{k{N}_{s,k}}+\frac{{\rho }^{\mathrm{‘}}}{k{N}_{s,k}{\sigma }_n^2}\left[\begin{array}{cc}{\mathsf{\Sigma }\mathbf{\text{'}}}_{kc,1}^2& 0\\ 0& {\mathsf{\Sigma }\mathbf{\text{'}}}_{kc,2}^2\end{array}\right]\left[\begin{array}{cc}\overline{Q}& 0\\ 0& 0\end{array}\right]\right|\right) \\ ={\log }_2\left(\left|{\mathit{I}}_{k{N}_{s,k}}+\frac{{\rho }^{\mathrm{‘}}}{k{N}_{s,k}{\sigma }_n^2}{\mathsf{\Sigma }\mathbf{\text{'}}}_{kc,1}^2\overline{Q}\right|\right) \\ ={\log }_2\left(\left|{\mathit{I}}_{k{N}_{s,k}}+\frac{{\rho }^{\mathrm{‘}}}{k{N}_{s,k}{\sigma }_n^2}{\mathsf{\Sigma }\mathbf{\text{'}}}_{kc,1}^2\overline{Q}\right|\right) \\ ={\log }_2\left(\left|{\mathit{I}}_{k{N}_{s,k}}+\frac{{\rho }^{\mathrm{‘}}}{k{N}_{s,k}{\sigma }_n^2}{\mathsf{\Sigma }\mathbf{\text{'}}}_{kc,1}^2\right|\right) \\ +{\log }_2\left(\left|{\mathit{I}}_{k{N}_{s,k}}-{\left({\mathit{I}}_{k{N}_{s,k}}+\frac{{\rho }^{\mathrm{‘}}}{k{N}_{s,k}{\sigma }_n^2}{\mathsf{\Sigma }\mathbf{\text{'}}}_{kc,1}^2\right)}^{-1}\frac{{\rho }^{\mathrm{‘}}}{k{N}_{s,k}{\sigma }_n^2}{\mathsf{\Sigma }\mathbf{\text{'}}}_{kc,1}^2\left({\mathit{I}}_{k{N}_{s,k}}-\overline{Q}\right)\right|\right) \\ \begin{array}{c}a\\ \approx \end{array}{\log }_2\left(\left|{\mathit{I}}_{k{N}_{s,k}}+\frac{{\rho }^{\mathrm{‘}}}{k{N}_{s,k}{\sigma }_n^2}{\mathsf{\Sigma }\mathbf{\text{'}}}_{kc,1}^2\right|\right)-tr\left({\mathit{I}}_{k{N}_{s,k}}-{\mathit{V}\mathbf{\text{'}}}_{kc,1}^H{\mathit{F}}_{RF}{\mathit{F}}_{BB,k}{{\mathit{F}}_{RF}}^H{{\mathit{F}}_{BB,k}}^H{\mathit{V}\mathbf{\text{'}}}_{kc,1}\right) \\ ={\log }_2\left(\left|{\mathit{I}}_{k{N}_{s,k}}+\frac{{\rho }^{\mathrm{‘}}}{k{N}_{s,k}{\sigma }_n^2}{\mathsf{\Sigma }\mathbf{\text{'}}}_{kc,1}^2\right|\right)-\left(k{N}_{s,k}-{‖{\mathit{V}\mathit{\text{'}}}_{kc,1}^H{\mathit{F}}_{RF}{\mathit{F}}_{BB,k}‖}_F^2\right) \end{gathered}$$

(75)

Where $\begin{array}{c}a\\ \approx \end{array}$ denotes the two approximations ${\left({\mathit{I}}_{k{N}_{s,k}}+\frac{{\rho }^{\mathrm{‘}}}{k{N}_{s,k}{\sigma }_n^2}{\mathsf{\Sigma }\mathbf{\text{'}}}_{kc,1}^2\right)}^{-1}\frac{{\rho }^{\mathrm{‘}}}{k{N}_{s,k}{\sigma }_n^2}{\mathsf{\Sigma }\mathbf{\text{'}}}_{kc,1}^2\approx {\mathit{I}}_{k{N}_{s,k}}$. When the $\ddot{Y}$ values are not very large, ${\mathrm{l}\mathrm{o}\mathrm{g}}_2\left|{\mathit{I}}_{k{N}_s}-\ddot{Y}\right|\approx {\mathrm{l}\mathrm{o}\mathrm{g}}_2\left(1-tr\left(\ddot{Y}\right)\right)\approx -tr\left(\ddot{Y}\right)$ , $\overline{Q}={\mathit{V}\mathit{\text{'}}}_{kc,1}^H{\mathit{F}}_{RF}{\mathit{F}}_{BB,k}{{\mathit{F}}_{RF}}^H{{\mathit{F}}_{BB,k}}^H{\mathit{V}\mathit{\text{'}}}_{kc,1}$ , and $\ddot{Y}={\mathrm{l}\mathrm{o}\mathrm{g}}_2\left|{\mathit{I}}_{k{N}_{s,k}}-\overline{Q}\right|$.

The system spectral efficiency is maximizing ${‖{\mathit{V}\mathit{\text{'}}}_{kc,1}^H{\mathit{F}}_{RF}{\mathit{F}}_{BB,k}‖}_F$ , which represents the optimal beamforming matrix ${\mathit{F}}_{opt}={\mathit{V}\mathit{\text{'}}}_{kc,1}$ and ${\mathit{F}}_{RF}{\mathit{F}}_{BB,k}$ the chord distance in the Grassmann manifold. The problem can be expressed as ${‖{\mathit{F}}_{opt}-{\mathit{F}}_{RF}{\mathit{F}}_{BB,k}‖}_F$. The beamforming optimization problem for a millimeter wave massive MIMO UAV is expressed as follows:

$$\left({\mathit{F}}_{RF}^{opt},{\mathit{F}}_{BB}^{opt}\right)=arg\underset{{\mathit{F}}_{\mathit{RF}}{,\mathit{F}}_{\mathit{BB}}}{\min }{‖{\mathit{F}}_{opt}-{\mathit{F}}_{RF}{\mathit{F}}_{BB,k}‖}_F$$

$$s_{.}{t}_{.}{\mathit{F}}_{RF}ϵ{\mathcal{F}}_{RF}$$

$${‖{\mathit{F}}_{RF}{\mathit{F}}_{BB,k}‖}_F^2={kN}_s$$

(76)

Because the constraint ${\mathit{F}}_{RF}ϵ{\mathcal{F}}_{RF}$ , that is${\mathit{F}}_{RF}$ the constant mode constraint on the column elements makes the problem non-convex and difficult to find an optimal solution. There is a need to improve the performance of the system while reducing the complexity and solving the non-convex optimization problem, so a hybrid beamforming optimization method based on PE-AltMin is introduced.

5.2. Digital and analogue beamforming matrix solutions

This section proposes an alternating minimization method to solve the above problem, and such methods have been successful in solving optimization problems on two variables. The optimization equation (76) is decoupled by keeping one variable fixed and optimizing iteratively the other variable. The optimization equation (76) is rewritten as the following equation:

$$\underset{{\mathit{F}}_{\mathit{BB}}}{\min }{‖{\mathit{F}}_{opt}-{\mathit{F}}_{RF}{\mathit{F}}_{BB,k}‖}_F$$

$$\underset{{\mathit{F}}_{\mathit{RF}}}{\min }{‖{\mathit{F}}_{opt}-{\mathit{F}}_{RF}{\mathit{F}}_{BB,k}‖}_F^2$$

$$s_{.}{t}_{.} \left|{\left({\mathit{F}}_{RF}\right)}_{i,j}\right|=1,\forall i,j$$

(77)

Two sets of alternating minimization methods (AltMin, Alternating Minimization) are defined in the literature [35]. One of them is a method for high complexity manifold optimization. In this case, the above problem can be solved iteratively by defining a Riemannian manifold. We first use the conjugate line-step gradient method to determine the analog beamformer under the constraints, after which we use the least squares method to find the digital beamformer. The first using the conjugate line step gradient method to determine the analogue beamformer under constraints, and the later using the least squares method to find the digital beamformer. The least-squares method is used to find the digital beamformer. The MO-AltMin (Manifold Optimization Alternating Minimization) method offers better performance as it optimizes both variables simultaneously. However, the complexity of the MO-AltMin method is relatively high in that the update of the simulated beamformer in each iteration involves a line search method, the conjugate gradient method. Therefore, the nested loops in the MO-AltMin method slow down the overall solution process. This requires the study of a hybrid beamforming method with low computational complexity and low performance loss. A low-complexity PE-AltMin is designed on the basis.

The constraint that the columns about the digital beamformer should be orthogonal is expressed as:

$${\mathit{F}}_{BB,k}^H{\mathit{F}}_{BB,k}=\alpha {\mathit{F}}_{DD,k}^H\alpha {\mathit{F}}_{DD,k}={\alpha }^2{\mathit{I}}_{{kN}_s}$$

(78)

Where ${\mathit{F}}_{DD,k}$ is a matrix with the same dimension as ${\mathit{F}}_{BB,k}$.Your matrix with the same dimension. The orthogonality constraint provides the possibility that the analogue beamformer can discard the product form with ${\mathit{F}}_{BB,k}$ product form, which will help to greatly simplify the design of the analogue beamformer. The objective function in (77) greatly simplifies the design of the analog beamformer when the AltMin method is applied. Since the matrix F_RF removes the product form of ${\mathit{F}}_{BB,k}$, the ${\mathit{F}}_{RF}$ closed form solution is expressed as:

$$arg\left({\mathit{F}}_{RF}\right)=arg\left({\mathit{F}}_{opt}{\mathit{F}}_{DD,k}^H\right)$$

(79)

Where $arg\left({\mathit{F}}_{RF}\right)$ generates a new matrix containing the phase of the${\mathit{F}}_{RF}$ matrix. The results show that the ${\mathit{F}}_{RF}$ the phase can be extracted from the equivalent beamformer ${\mathit{F}}_{opt}{\mathit{F}}_{DD,k}^H$. This closed form of the solution can also be seen as ${\mathit{F}}_{opt}{\mathit{F}}_{DD,k}^H$ on the feasible set of the simulated beamformer ${\mathcal{F}}_{RF}$ on the Euclidean projection. However, this does not affect the complexity of the digital beamformer and modifies equation (77) to (80):

(80)

Based on the low-complexity PE-AltMin method, we reduce the search range of the digital beamformer. Equation (77) is expressed as:

$${\mathit{F}}_{DD,k}{\mathit{F}}_{RF}={\mathit{F}}_{opt}$$

(81)

In equation (80), since the combined optimization problem eliminates the need for ${\mathit{F}}_{DD,k}$ the other solutions, the complexity of the PE-AltMin method can be better reduced. Furthermore, instead of randomly selecting the phase during the initialisation of the PE-AltMin method, equation (79) is used to select ${\mathit{F}}_{RF}$. Finally, the final solution of the digital beamformer can be obtained by the least squares correction given in Equation (82).

$${\mathit{F}}_{DD,k}={\left({\mathit{F}}_{RF}^H{\mathit{F}}_{RF}\right)}^{-1}{\mathit{F}}_{RF}^H{\mathit{F}}_{opt}$$

(82)

In millimeter-wave massive MIMO of UAV, the number of antennas at the transmitting would be large, where ${\mathit{F}}_{DD,k}ϵ{\mathbb{C}}^{N_{RF}^T\times N_{S,k}}$ ,${\mathit{F}}_{RF}ϵ{\mathbb{C}}^{N_T\times N_{RF}^T}$ , ${\mathit{F}}_{opt}ϵ{\mathbb{C}}^{N_T\times N_{S,k}}$. The time consumed and the corresponding complexity of the matrix ${\mathit{F}}_{opt}{\mathit{F}}_{RF}$ performing SVD gradually increases. So the projection approximation subspace tracking method with lower complexity and better real-time is designed to estimate only ${\mathit{F}}_{opt}{\mathit{F}}_{RF}$ of the right singular matrix $\overbrace{S}$. This avoids the need for SVD to obtain the entire left and right singular matrices and the computation of the singular value matrix.

In the projection approximation subspace tracking method, the ${\mathit{F}}_{opt}{\mathit{F}}_{RF}$ as the data sample. ${\mathit{F}}_{opt}{\mathit{F}}_{RF}$ is selected a row in the data sample vector $\overbrace{x}$ . When inputting the $\overbrace{j}\left(1\le \overbrace{j}\le N_T\right)$ row ${\overbrace{x}}_{\overbrace{j}}$ is entered, it is possible to obtain ${\mathit{F}}_{opt}{\mathit{F}}_{RF}$ of $\overbrace{s}$ (1$\le \overbrace{s}\le \overbrace{S}$ ) column${\overbrace{\mathit{B}}}_j\left(:,\overbrace{s}\right)$. ${\overbrace{\mathit{B}}}_j\left(:,\overbrace{s}\right)$ can be expressed as the following equation:

$${\overbrace{\mathit{B}}}_{\overbrace{j}}\left(:,\overbrace{s}\right)={\overbrace{\mathit{B}}}_{\overbrace{j}-1}\left(:,\overbrace{s}\right)+\frac{1}{d_{\overbrace{j}}\left(\overbrace{s}\right)}\left[x_{\overbrace{j}}-{\overbrace{\mathit{B}}}_{\overbrace{j}-1}\left(:,\overbrace{s}\right)\overbrace{y}\right]{\overbrace{y}}^{\mathrm{*}}$$

(83)

Where $\overbrace{y}={\overbrace{\mathit{B}}}_{\overbrace{j}-1}^H\left(:,\overbrace{s}\right){\overbrace{x}}_{\overbrace{j}}$ is the intermediate variable, and ${\overbrace{y}}^{\mathrm{*}}$ denotes the $\overbrace{y}$ the conjugate. ${d\mathrm{\text{'}}}_{\overbrace{j}}\left(\overbrace{s}\right)={d\mathrm{\text{'}}}_{\overbrace{j}-1}\left(\overbrace{s}\right)+{\left|{\overbrace{y}}_{\overbrace{j}}\right|}^2$ is the step size used for the correction. For the method steps can be described as follows:

Algorithm 1: Projection Approximation Subspace Algorithm

Applying the above method to the improved hybrid beamforming optimization method of finding PE-AltMin for ${\mathit{F}}_{opt}^H{\mathit{F}}_{RF}^{\left(\stackrel{`}{k}\right)}$ the SVD calculation, then the specific method based on PE-AltMin can be expressed as follows:

Algorithm 2: Improved hybrid beamforming optimization algorithm based on PE-AltMin

6. Simulation results and analysis

For evaluation, we perform Monte-Carlo simulations. We present numerical results forOutage probability, the transmission rate and BER (Bit Error Rate) performance of FSO links for various numbers of transmit/receive apertures and correlation values. We consider a FSO system with a receive aperture of size $D_0$ = 5cm and the wavelength of $\lambda$ = 1.55μm. The link distance is assumed to be $L$ = 2km. The correlation length can be therefore approximated as $d_0$ ≈ $\sqrt{\lambda L}$ = 5.5cm. To validate the performance of these models millimete-wave massive MIMO/FSO in UAV system, the carrier frequency is 78GHz. The grid size of the dictionary matrix is set to $G_R=2N_R=64$,$G_T=2N_T=128$. The maximum number of channel paths is set to $Q_{max}=5$. The UAV speed is 0-80 km/h.

We compared the performance of four methods in terms of FSO transmission rate, as shown in Figure 5. As can be seen, the dual-hop FSO in UAV system with BiGRU-Attention neural network method offers the best performance in terms of the average FSO transmission rate. For instance, at SNR=20dB, using BiGRU-Attention network results in the average FSO transmission rate of 1050Mbps, 900Mbps, 700Mbps, 600Mbps, compared with the DHPL-ML, LDPC-MIMO/FSO, and RA-FSO-RF, respectively. The results indicate that the application of BiGRU-Attention neural network model in FSO system can effectively solve the optical PE problem caused by UAV jitter. Although other methods can improve optical communication performance, the ability of FSO PE to handle complex jitter in FSO environments is limited.

We compared the BER of four methods. As shown in Figure 6, the dual-hop FSO in UAV system with BiGRU-Attention neural network method offers the best performance in terms of the average BER. When the SNR=25dB, using BiGRU-Attention network results in a decrease in the average BER of ${10}^{-5}$, ${10}^{-5}$, ${10}^{-5}$, ${10}^{-5}$, compared with the DHPL-ML, LDPC-MIMO/FSO, and RA-FSO-RF, respectively. Obviously, the BiGRU-Attention network model achieves lower BER compared to other methods. In addition, all considered methods exhibit a trend of decreasing error rate with increasing SNR. This also reveals the significant impact of PE of UAV in these systems.

We investigate the performance of the system under consideration based on the probability of interruption. Figure 7 shows the relationship between interruption probability and SNR ratio under different methods. Under the same SNR ratio conditions, the BiGRU-Attention network performs the best. For instance, at SNR=20dB, using BiGRU-Attention network results in outage probability of ${10}^{-4}$, ${10}^{-4}$, ${10}^{-4}$, ${10}^{-1}$, compared with the DHPL-ML, LDPC-MIMO/FSO, and RA-FSO-RF, respectively.

According to the simulation results shown in Figure 10, the spectral efficiency of the proposed method in this paper reaches 14bps/Hz with SNR=20 dB, which is the best performance. In comparison, the spectral efficiency of the RBOP method is 12.5 bps/Hz, the REHS method is 11 bps/Hz, the ZF-HBF method is 8 bps/Hz, and the Hy-BD method has a spectral efficiency of only 6 bps/Hz, which is the worst performer.

We compared the BER of five methods, as shown in Figure 11. When the SNR=25dB, the bit error rate of the proposed method is ${10}^{-5}$, which is better than other methods. The bit error rate of the RBOP method is ${0.5\times 10}^{-4}$, the bit error rate of the REHS method is ${0.7\times 10}^{-4}$, the bit error rate of the ZF-HBF method is ${1.2\times 10}^{-4}$, and the bit error rate of the Hy-BD method is the highest, reaching ${5\times 10}^{-4}$.

Figure 12 shows that when the maximum Doppler shift is 8kHz, the BER of the proposed method reaches ${10}^{-3}$, which is the best performance. The BER of the RBOP method is $5.5\times {10}^{-3}$, because the RBOP method has high channel requirements and cannot effectively suppress multipath interference, which affects the improvement of system performance. The REHS method has BER of ${10}^{-2}$. The REHS method requires a lot of computation and processing, which reduces the performance of the system. The BER of the ZF-HBF method is $5.5\times {10}^{-2}$. The Hy-BD method has the highest BER of ${10}^{-1}$. The Hy-BD method performed the worst.

7. Summary

In this paper, a Tensor-train decomposition based hybrid beamforming for millimeter-wave massive MIMO/FSO in UAV with RIS Networks is investigated. Firstly, the system channel model of millimeter-wave massive MIMO/FSO in UAV with RIS is established. The application of BIGRU (Bidirectional Gated Recurrent Unit)-Attention neural net-work model in an FSO system can effectively solve the problem of light PE caused by UAV jitter. The multidimensional received signal of millimeter-wave massive MIMO in UAV with RIS is represented as a low-rank tensor signal by Tensor-train decomposition. The problem of maximizing the system spectral efficiency is decomposed into two subproblems of optimizing the RIS phase shift matrix and solving the hybrid beamforming matrix. The hybrid channel matrix is decomposed by SVD. The RIS phase shift matrix is optimized. Finally, the hybrid beamforming matrix is solved by decoupling the transceiver and transmitter using an improved hybrid beamforming method based on PE-AltMin.