Low complexity channel estimation algorithm using paired spatial signatures for UAV 3D MIMO systems

This paper proposes a low complexity channel estimation algorithm for unmanned aerial vehicle three dimension multi-user multiple-input-multiple-output (3D MIMO) systems with the uniform planar array (UPA) at base station using paired spatial signatures. With the aid of antenna array theory and array signal processing, 3D channel is firstly modeled based on the angles between the direction of arrival along x- and y-axis of the UPA. And 3D MIMO channels can be projected onto the x- and y-directions, respectively. Then, channel estimation for multi-user uplinks using small amount of training resources is divided into two phases. At the first uplink preamble phase, each user is assigned the orthogonal pilot, and the paired spatial signatures and optimal rotation angle of each user through the same pilot sequence are obtained. We also propose a user grouping strategy based on three-dimension angle-division multiple access (3D-ADMA) to ensure that the user's spatial signatures do not overlap. At the second phase during several coherence times, the same pilot sequence within a group and orthogonal pilot sequences between groups are assigned, then, the channel state information of the user's x- and y-directions are recovered by the paired space signatures and optimal rotation angle of each user obtained in the preamble phase, respectively. And dynamically updating the user's paired spatial signatures and optimal rotation angle utilizes the obtained channel parameter of x- and y-directions. Finally, the channel parameter of the x- and y-directions are reconstructed by the updated user's space signatures and the optimal rotation angle, and the 3D MIMO channel estimation is obtained through the Kronecker product. Compared with the conventional channel estimation method of the 3D MIMO system under UPA using a low-rank model, the proposed methods reduce the computational complexity without degrading the estimated performance to a large extent. Furthermore, it is carried out with limited training resources, and the pilot resource overhead of the system is greatly reduced by the 3D-ADMA packet and the two-stage pilot allocation. Simulation results verified that the proposed algorithm is effective and feasible.


Introduction
During recent years, the research on 5th-generation (5G) technology has developed at the speed of blowout.Massive multiple-input-multiple-output (Massive MIMO) system [1] have drawn considerable interest from academia and industry.Massive MIMO is a new technique that employs hundreds or even thousands of antennas at many types of base stations (BSs), such as unmanned aerial vehicle (UAV)-based BS, to simultaneously serve multiple single-antenna users, and has been widely investigated for its numerous merits, such as high spectrum and energy efficiency, high spatial resolution, and reducing network interference [1][2][3].Considering the arrangement of the antennas, the required antenna panel will be large [4] if the large-scale antenna array is only deployed in the horizontal dimension.Consequently, placing the antenna in a two-dimensional grid can reduce the size of the antenna panel, which is called three-dimension multipleinput multiple-output (3D-MIMO).3D-MIMO can not only utilize the spatial freedom of large-scale transmit antennas but also adjust the direction of the transmit beam in horizontal and vertical dimensions, which improves spatial resolution, improves signal power, and reduces inter-cell interference [3].UAV communications have been widely used in wireless communications recently such as communication relay, information dissemination, data collection and so on, due to the controllable mobility and everdecreasing manufacturing cost [5,6].Sharma et al. [7] concentrated on a simple path loss and shadow fading channel model that is commonly used to describe the propagation between an aerial base station and a user on the ground.Jiang et al. [8] proposed a 3-D elliptic-cylinder unmanned aerial vehicle (UAV) multiple-input-multiple-output (MIMO) channel model for air-to-ground communication environment, which expressed by the space correlation functions, provides a novel and practical approach to investigate UAV-MIMO channels and design vehicular communication systems.Sudheesh et al. [9] considered the application of IA in a high altitude platform (HAP) to ground station (GS) communication, and the application of IA is proposed for a generalized channel in a HAP-to-GS communication link that takes into account angle-ofdeparture and angle-of-arrival at the transmitter and at the receiver, respectively.
In a variety of communication application scenarios, the accurate channel state information (CSI) between BS and users is very important.The CSI acquisition has been recognized as very challenging task for 3D MIMO systems, due to the high dimensionality of channel matrices as well as the resultant uplink pilot contamination, prohibitive computational complexity and so on [10].To solve these problems, some researchers have used the sparsity of the Massive MIMO channel to estimate the channel with a small number of observations, thereby reducing computational complexity and pilot overhead [11][12][13][14].The author of [11] proposed the Massive MIMO channel estimation problem as a compressed sensing (CS) problem.However, it is difficult to obtain the sparseness of the channel in actual operation.It is generally assumed that the channel sparse energy level is known, i.e., the number of non-zero elements of the channel impulse response (CIR) is known [12,13].The authors of the literature [14] reconstructed the sparse channel of the Massive MIMO system with a small amount of observation data through the Bayesian learning method.In addition to the above research directions, some methods for channel estimation using the low rank of the channel covariance matrix are proposed [15][16][17].In the literature [15][16][17], it is assumed that each user arrives at a narrow angle, and then the eigenvalue decomposition of the channel covariance matrix can effectively reduce the channel dimensions [18].
The author of [19] proposed a discrete-Fourier-transform-based (DFT-based) spatial basis expansion model (SBEM) method in multi-user Massive MIMO systems.The method transforms the uplink/downlink channel into a spatial domain capable of exhibiting sparsity by using a DFT transform matrix according to the physical characteristics of the uniform linear array.And avoiding the problem of uplink channel pilot interference by grouping and scheduling the user's spatial signatures.And grouping and scheduling are based on angle division multiple access (ADMA) [20].In addition, the channel reciprocity of the TDD system is used to simplify the downlink training sequence through the uplink spatial information.The SBEM method does not require channel statistics compared with the previously described channel estimation method using channel sparsity.Xie et al. [21] proposed to reconstruct the channel covariance matrix by estimating the angle information and power angle spectrum of the channel based on the SBEM method, and finally enhance the accuracy of channel estimation by the minimum mean square error (MMSE) estimator.And along with the characteristics of a uniform linear array, there is an optimal rotation angle that can make the receive power of the user become more concentrated after the channel vector is DFT-transformed.Thus, the channel information of the user becomes sparser, and the channel information can be estimated by less observation information.
Author in [22] proposed a channel estimation algorithm based on 2D discrete Fourier transform (2D-DFT) for indoor 60 GHz massive MIMO systems via array signal processing.The algorithm requires a small amount of training overhead, which greatly reduces the training overhead and feedback cost, and increases the number of user terminals that the systems can serve by utilizing the different spatial information.Fan et al. [23] proposed a channel estimation algorithm based on 2D-DFT for millimeter-wave Massive MIMO system via the physical characteristics of the antenna array.Firstly, the arrival angle information of each user's different paths is estimated, then the accuracy of the estimation is enhanced by the angle rotation technique, and the gain information of each path is estimated.Finally, the channel is reconstructed.However, as the number of antennas increases in a uniform planar array, the required computational complexity will increase exponentially, which undoubtedly increases the computational overhead of the system.
In the channel estimation problems of Massive MIMO systems, angle information plays a very important role.Therefore, accurate angle information estimation methods are urgently needed.In the field of array signal processing, angle estimation algorithms based on traditional high-resolution subspaces, such as multiple signal classification (MUSIC), estimating sign parameters via rotational invariant technique (ESPRIT) and so on, are of great interest due to their high-resolution angle estimation algorithms [24][25][26].These methods have been extensively studied in traditional Massive MIMO systems and 3D MIMO systems [27][28][29][30].Both the MUSIC algorithm and the ESPRIT algorithm utilize the covariance matrix of the received signal and perform singular value decomposition (SVD) operations on the covariance matrix.The computational complexity of SVD operations under large-scale antenna arrays can be high and the category of blind estimation does not fully utilize a priori information in wireless communications.
In the field of antenna array signal processing, the authors in [31] proposed an improved model to improve direction-of-arrival estimation accuracy, not only for uniform planar array (UPA) but also for sparse L-shaped array.The author converts the twodimensional UPA signal model represented by the signal elevation angle and azimuth angle into the two-dimensional UPA signal model represented by the electrical angles in x-and y-directions of the signal.This simplifies the UPA signal model and is easier and less computationally complex for 2D direction-of-arrival problems.
Inspired by the above literature, this paper focuses on the research of UAV 3D MIMO systems.Simplify the channel model and project the received signal onto the uniform linear array (ULA) on the x-and y-axis, respectively.By simplifying the channel model, the data projected to the ULA is only related to one angle of the 3D space, then DFT and angle rotation techniques are used to achieve super-resolution estimation, and the calculation complexity is greatly reduced.Finally, the 3D channel is obtained through the Kronecker product.Through the grouping and pilot allocation strategy, the orthogonal pilot overhead is reduced, and the problem of pilot interference caused by non-orthogonal pilots is solved.
Notations Small and upper bold-face letters refers to column vectors and matrices, respectively; the superscripts (•) H , (•) T , (•) , (•) −1 stand for the conjugate-transpose, trans- pose, conjugate, inverse of a matrix, respectively; [A] i,j is the ( i, j th entry of A ; diag{a} denotes a diagonal matrix with the diagonal element constructed from a; ≜ represents new definition; − → (A) denotes the vectorization of A ; [h] B , indicates the sub-vector of h by keeping the elements indexed by B ; [H] B , stands for the sub-matrix of H by collecting the rows indexed by B ; [H ] B denotes the sub-matrix of H by collecting the columns indexed B.

Method
This paper considers the UAV 3D MIMO scenario.A UPA antennas of M N is configured at the BS to serve all users configuring a single antenna.Simplify the 3D MIMO channel model using the idea of [31], channel modeling by the azimuth and elevation angles is transformed into that channel modeling uniquely determined by the electrical angles in x-and y-directions of the signal.Under the simplified channel model, the 3D channel can be projected onto the x-and y-directions, respectively.The angle information of the channel projected onto the x-direction is only related to the angle of the x-direction, and the angle information of the channel projected onto the y-direction is only related to the angle of the y-direction.Therefore, the information projected onto the x-direction and the y-direction is estimated separately, and finally the 3D channel is obtained through the Kronecker product.Since the x-direction and the y-direction are composed of M and N antennas, respectively, the initial DOAs can be estimated from the physical properties of ULA through the DFT matrices of the M and N points, respectively.The angle rotation technique is then used to further enhance the estimation precision.
Compared with the 2D-DFT methods, the proposed method only requires twice DFT matrix operations, which greatly reduces the computational complexity.Inspired by the idea of [20], we propose a 3D-ADMA grouping and scheduling strategy in 3D space, so that the angle between the x-direction and the y-direction does not overlap and the certain guard interval is satisfied.And the pilot interference caused by the use of nonorthogonal pilot sequences within the cell is reduced.It is worth noting that there is a problem of user angle paired when acquiring the initial angle of arrival of the user.At this stage, the angles obtained in the x and y directions are paired by the same pilot sequence, thereby obtaining the user's paired spatial signatures and optimal rotation angles.It should be noted that the user's paired spatial signatures refer to an index set of DFT points of channel energy concentration after DFT and rotation operations on the user channel vector.And each user has two index sets.Then all users are grouped by 3D-ADMA, the same pilot sequence is assigned within the same group, and orthogonal pilot sequences are used in different groups.The two-stage pilot allocation scheme can greatly reduce pilot overhead.

Contributions
The main contributions of the paper are the following: 1. Utilizing the method of signal modeling in reference [31], we simplified the channel modeling of UAV 3D MIMO system based on the angles between incident signal and x-axis and y-axis of UPA.In this paper, the received signal is projected to x-axis and y-axis respectively, the angles along x-axis and y-axis can be obtained respectively, and 3D MIMO channel is generated by Kronecker product.2. In reference [18][19][20][21], the channel estimation and data transmission of Massive MIMO system are well studied.In this paper, the ADMA user grouping strategy in reference [20] is introduced into UAV 3D MIMO system, and a 3D-ADMA user grouping strategy is proposed to save pilot overhead.The idea of SBEM channel estimation method based on DFT in reference [18,21] is used to update the angle information of users and reconstruct 3D MIMO channel.In order to solve the problem of matching the angles along x-axis and y-axis, the paper uses the same orthogonal pilot to pair the projected data to obtain the spatial characteristics and the optimal rotation angle.This paper proposes a new channel estimation method, which includes modeling, grouping, using pairwise spatial features and optimal rotation angle and SBEM channel estimation scheme, will generate fresh insight into channel estimation of 3D-MIMO system.
The rest of the paper is organized as follows.In Sect.3, the system model of 3D MIMO system and channel model are described.In Sect.4, a two-stage channel estimation algorithm based on twice DFT is introduced.Numerical results and discussions are then provided in Sect. 5. Finally, conclusions are drawn in Sect.6.

System models
Let us consider such a 3D MIMO scenario, where BS equipped with antennas serves K single-antenna users.Nowadays, most of BS are built on high buildings or very high signal towers in the macrocells.In this case, the signal sent by the user to the BS will be limited to a narrow range of incident angles.And let's enable UAV-based BS systems when traditional BS fails, which is assembled in a UPA.Therefore, let us consider the Saleh-Valenzuela (SV) spatial channel model, the antenna array uses a UPA, and the UPA signal model is shown in Fig. 1. M and N antennas are placed on the x-and y-directions, respectively, and the antenna spacing of the x-axis and the y-axis is d.

Transmitter model
In the uplink transmission stage, the received signal at the BS can be expressed as where Y ∈ C (M×N )×L is signal matrix received at the BS, L is pilot length sent by the user, h k ∈ C (M×N )×1 is channel vector between the BS and the kth mobile user, H ∈ C (M×N )×K is channel matrix between the BS and all users, S ∈ C K ×L is user's trans- mit pilot matrix, and N ∈ C (M×N )×L is the complex Gaussian noise matrix.

Channel model
Suppose the array antenna receives K user data, and let us define φ k ∈ [−π/2, π/2) and θ k ∈ [−π , π) as the signal elevation angle and azimuth of the kth user shown in Fig. 1, and the angle between the incident signal and the x-and y-directions are defined as α k and β k , respectively.From the geometric relationship, we can know that the relationship between these four angles can be expressed as

Then we can get
It can be seen from Eqs. ( 2) to ( 5) that the 3D MIMO channel can be uniquely determined not only by the signal elevation angle φ k and azimuth θ k , but also by the angles α k and β k between the incident signal and the x-and y-directions.This conversion pro- vides convenience for 3D MIMO channel modeling and signal processing.Therefore, 3D MIMO channel modeling can be expressed as where T is array manifold vector (AMV) or array T is AMV of ULA on the y-direction for β k,r , denotes the signal carrier wavelength, d is the antenna spacing and is assumed an uniform linear array with half-wavelength spacing (i.e., d = /2 ), z k,r denotes complex gain, and R is R rays in a narrow angular range.It should be noted that the vector h k in Eq. ( 1) is rearranged from matrix H k in Eq. ( 6), that is, h k vec(H k ).

Low complexity channel estimation algorithm
Figure 3 shows the transmission data phase diagram.First, the data transmission begins with the uplink preamble phase, and the initial paired spatial signatures and the optimal rotation angle are obtained by extracting the data received by the receiving end x-and y-directions antenna end (i.e., data projected to the x-and y-axis directions), respectively. ( The initial user grouping is then done by the initial paired spatial signatures.The transmission in the subsequent coherence time is used for channel estimation of the uplink through limited pilot resources.The operation here is also to extract the data of the receiving end portion for channel estimation processing, and dynamically update the angle information.
At the same time, it is detected whether the user angle exceeds the system threshold.When the angle information exceeds the threshold, the uplink preamble phase is re-executed and the users are grouped and scheduled.

Uplink preamble phase
In the uplink preamble phase, all K users obtain their initial CSI by the conventional Least squares (LS) method.Prior to this, the data received by the antennas on the x-and y-directions are first extracted at the receiving end of the BS Y x and Y y , respectively, where The initial CSI of each user's x-and y-directions is estimated using the LS method, denoted as ĥini k,x and ĥini k,y , k = 1, 2, . . ., K .The initial angle state information corresponding to the user's x-and y-directions is then extracted from the CSI.
For the sake of simplicity, let us assume that each user arrives at the BS with only one ray and the gain is constant at one.The receiving matrices projected onto the x-and y-directions are denoted as where and y-directions uplink channel matrix for all K users, respectively.

S
s T 1 , s T 2 , . . ., s T K T ∈ C K ×L denotes orthogonal pilot sequence matrix for all K users, and N x , N y are Gaussian noise matrix with CN (0, 1) elements.Then the channel vectors h k,x and h k,y of the x-and y-directions can be estimated by the LS method, respectively, and be expressed as (7) where n k,x and n k,y are Gaussian noise vectors of x-and y-directions.In the case of multi- users, the pilots allocated by each user are orthogonal in the preamble phase.Under such conditions, s k s H k = 1 , the application of the LS method can be expressed as formula (9) and formula (10).Repeating the ( 9) and (10) operations can obtain initial channel estimates for the x-and y-directions for all users.
The next step is to obtain initial angle information for each user's x-and y-directions through DFT and angle rotation techniques.We first define two normalized DFT matrix M pq , p, q = 0, 1, . . .M − 1 and Define: where ψ k,x = diag 1, e jψ k,x , . . ., e j(M−1)ψ k,x , ψ k,x ∈ [−(π/M), π/M] and ψ k,y = diag 1, e jψ k,y , . . ., e j(N −1)ψ k,y , ψ k,y ∈ [−(π/N ), π/N ] are spatial rotation parameter projected onto the x-and y-directions, respectively [19].When the optimal spatial rotation angle is found, the channel power can be more concentrated on a small number of DFT points, thereby reducing energy leakage and obtaining more accurate DOA estimation.When the channel has only one path composed of one ray, hro k,x and hro k,y have only one non-zero element p 0 and q 0 after the rotation operation, and as B ro k,x = {p 0 } and B ro k,y = {q 0 } are the paired spatial signatures of the kth user projected on the x-and y-directions, respectively.However, B ro k,x and B ro k,y in a multi-ray composition path are multi-element index set.Therefore, the kth user paired spatial signatures B ro k,x , B ro k,y and the optimal rotation angle ψ k,x , ψ k,y can be obtained from ĥini k,x and ĥini k,y by DFT and angular rotation techniques, respectively.Initial angle information αini k,x and βini k,x in which the kth user matches the x-and y-directions can then be obtained by the Eqs.( 13) and ( 14), respectively.
It should be noted that it may be expensive to allocate one orthogonal pilot sequence pilot overhead for each user in the uplink preamble phase in the whole system, but the operation of this phase is not used frequently, so the pilot overhead of the whole system is relatively small.The operation of the uplink preamble phase will only be restarted when the user's angle changes greatly.The users are grouped and scheduled.(11)

3D space grouping strategy
All users are grouped by the paired spatial signatures B ro k,x , B ro k,y (or initial angle information αini k,x , βini k,x ) obtained in the uplink preamble phase.The spatial signatures in the x-and y-direction are not overlapped, respectively, and all maintain a certain guard interval, and users who meet these two conditions are divided into the same group, called 3D-ADMA.Divide K users into G groups, denoted {[U 1 , U 2 , . . ., U G ]} .Grouping must satisfy (15)  expression.
where define y , The setting of the thresholds x and y depends on the tolerance of the multiplexed pilots in the system [21].

Intra-group channel estimation
The uplink preamble phase is a very short phase and may have only one coherence time or less than one coherence time.After grouping, the users are divided into G groups.In each subsequent coherence time, each group in the uplink channel transmission shares one pilot sequence, called the Intra-group pilot multiplexing, and the orthogonal pilot is used between the groups.This will avoid inter-group interference and greatly reduce pilot overhead.Intra-group pilot multiplexing will cause pilot interference, but the paired spatial signatures of the users in the group are not overlapping, so that the users in the group implement orthogonal transmission.Since the channel coherence time is in the millisecond level and the building around each user does not physically change in a short time, the user's DOA can be seen as having not changed in dozens of coherent times.For example, if coherence time equals 5 ms, the maximum displacement for a user moving at 80 km/h in ten coherence times will be only 1.11 m, which means that the resultant angle variation seen by the BS is negligible [21].Therefore, the uplink channel estimation after the uplink preamble phase can use the paired spatial signatures B ro k,x and B ro k,y obtained in previous coherent time (or in uplink preamble phase).
For the convenience of analysis, we will analyze all users in the U 1 group as an example.Orthogonal pilot sequences assigned to the U 1 group by s 1 , and satisfy s 1 s H 1 = ρ u , where ρ u is the total uplink training signal to noise ratio (SNR).Then, in the nth coherent time after the end of the uplink preamble phase, the received signal of the BS from all users in the U 1 group can be expressed as where N U 1 is Complex Gaussian white noise matrix with CN (0, 1) elements.Orthogonal pilot sequences are used between groups, so there is no inter-group interference during (15 transmission.Then, the received signals of the BS UPA in the x-and y-directions are extracted from the Eq. ( 16), respectively, as where N U 1 ,x and N U 1 ,y are Complex Gaussian white noise matrix.Then the uplink chan- nel estimates on the x and y-directions of all users in the U 1 group are represented by the conventional LS estimation algorithm as is the normal- ized noise vector in the x-and y-directions, respectively.Since the DOA of all users in the U 1 group do not overlap, the channel information of the kth user in the U 1 group can be extracted from the aliased h U 1 ,x (n) and h U 1 ,y (n) by the paired spatial signatures obtained in previous coherent time (or in uplink preamble phase).First do DFT and angle rotation techniques for h U 1 ,x (n) and h U 1 ,y (n) , and then extract the kth user's chan- nel information by using the paired spatial signatures B ro k,x (n − 1) , B ro k,y (n − 1) and opti- mal rotation angle ψ k,x (n − 1) , ψ k,y (n − 1) of the (n-1)th coherent time of the kth user.The DFT and angular rotation of the channel projected to the x-and y-directions of the nth coherence time can be expressed as [24] (17 Since B ro k,x (n − 1) and B ro i,x (n − 1) are not overlapping, and B ro k,y (n − 1) and B ro i,y (n − 1) are also non-overlapping, we can eliminate the second term after the first equal sign in Eqs. ( 23) and (24), respectively.Note: The spatial characteristics and the optimal rotation angle information of the (n − 1)th coherence time are used in calculating the channel estimation phase of the nth coherence time.This is because the DOA of each user in a short period of time changes very weakly.Then we can recover the channel information of the kth user projected to the x-and y-directions from ( 21) and ( 22), respectively, as Next, we recover the 3D MIMO channel of the kth user's nth coherence time through the Kronecker product.
When obtaining the CSI of the nth coherent time of the k user, we use spatial signatures and the optimal rotation angle of the previous coherence time.However, spatial signatures and the optimal rotation angle using the nth coherence time is more accurate for channel estimation than spatial signatures and the optimal rotation angle in the previous coherence time.That is, Therefore the paired spatial signatures B ro k,x (n) , B ro k,y (n) and the optimal rotation angle ψ k,x (n) , ψ k,y (n) of the nth coherent time update of the k user are obtained from ) ,respectively.Then the angle information αk,x (n) and βk,y (n) of the x-and y-directions projected can be estimated, which can be expressed as (22) k,y (n−1),: Then, the estimated angle information is used to determine whether the angle of each two users exceeds a threshold set by the system.If it is exceeded, the uplink preamble phase is re-executed, and user grouping and scheduling are resumed.
Next, we reconstruct the channel by the angles estimation αk,x (n) and βk,y (n) by updating the paired spatial signatures and the optimal rotation angle.For the convenience of analysis, the gain is assumed to be 1, so that the channel reconstructing the x-and y-directions can be expressed as Therefore, the kth user's nth coherent time 3D MIMO channel estimation is reconstructed by Kronecker product, expressed as Repeating the above operations in the U 1 group can reconstruct the channel informa- tion of all users in the U 1 group.By repeating the intra-U 1 operation in each group in the cell, we can complete the uplink channel estimation for all users in the cell.

Computational complexity
Compare the computational complexity of algorithm in the section.The computational complexity of the proposed method, the channels projected onto the x-and y-directions are processed separately using the SBEM method in [15] and generating a 3D MIMO channel by Kronecker product, method of 2D-DFT and angular rotation technique in [20] are expressed in terms of the required number of floating point operations (FLOPs) and shown in Table 1 and Fig. 4. Therein, K represents the number of (26) and

6MN
SBEM in [21] UL preamble LS method to Y DFT and spatial rotation method single antenna users served by the BS, G represents the number of groups, N g represents the number of users in each group, N ro represents to find the optimal number of rotation angles with angle rotation technology, N co represents the number of coherence times, N b represents the number of the user's spatial signatures, and one complex-valued mul- tiplication and one complex-valued addition require 6 and 2 FLOPs, respectively [32].As can be observed from Table 1 and Fig. 4, the complexity of method of 2D-DFT and angular rotation technique in [23] increased cubically with number of antennas of the BS.While the computational complexity of the proposed method and the channels projected onto the x-and y-directions are processed separately using the SBEM method in [18] increase only linearly with number of antennas.Since number of antennas is considerably larger than the other parameters in 3D MIMO systems and is much larger than the traditional Massive MIMO systems, method of 2D-DFT and angular rotation technique in [20] is the most computational complexity, and the computational complexity of the proposed estimation method is slightly higher than that of the method,which projected received signal to the x-and y-directions and processed separately using the SBEM method in [18], but greatly lower than method of 2D-DFT and angular rotation technique in [20].

Numerical results and discussions
In this section, we demonstrate the effectiveness of the proposed method through numerical example.we consider a multi-user 3D MIMO system with BS configuration UPA, and system parameters are chosen as M = 100 and N = 100, which are number of antennas in the x-and y-directions of the UPA, respectively, antenna spacing d = λ/2.K = 16 users are evenly distributed in the base station service range, and users are divided into 4 groups according to the user's distribution, and there are 4 users in (  6 Uplink MSE of proposed method as a function of SNR.Compares the performance of the proposed method with the channels projected onto the x-and y-directions are processed separately using the SBEM method in [18] then generating a 3D MIMO channel by Kronecker product, and method of 2D-DFT and angular rotation technique in [23] Figure 5 illustrates the MSE performance as a function of SNR for the proposed methods with lengths of pilot sequences L = 16 , L = 32 , and L = 64 , respectively.All users K = 16 are divided into 4 groups, and each group has 4 users.The system assigns a pilot sequence to group and guarantees the allocation of orthogonal pilot sequences between groups.The optimal number of rotation angles is N ro = 60 .It is shown in Fig. 4 that, as the length L of the pilot sequence increases, the performance of the uplink MSE also increases.And there is same error floor for all values of L. This phenomenon is not unexpected due to the truncation error of SBEM in [18] for the real channel and can be also observed in temporal Basis Expansion Model (BEM) [33].But, since the truncation error is only related with the effective expansion number (i.e., number of B k indexes), the error floors will remain the same for different L.
Figure 6 compares the performance of the proposed method with the channels projected onto the x-and y-directions processed separately using the SBEM method in [18] then generating a 3D MIMO channel by Kronecker product, and method of 2D-DFT and angular rotation technique in [23].K = 16 users are evenly distributed in the base station service range, and users are divided into 4 groups according to the user's distribution, and there are 4 users in each group, and L = 16 , N ro = 60 .It can be seen from the Fig. 5 that the estimated MSE performance of the proposed method is better than projected onto the x-and y-directions processed separately using the SBEM method in [18].And it is between the SBEM method in [18] then generating a 3D MIMO channel by Kronecker product and the method of 2D-DFT and angular rotation technique in [23].However, it can be seen from Fig. 4 that the computational complexity of the method of 2D-DFT and angular rotation technique in [23] is much larger than that proposed in this paper.Therefore, the method proposed in this paper greatly reduces the computing resource overhead without reducing the estimation performance.
Figure 7 illustrates the variation of angle MSE as SNR for the proposed method with different search times of angle rotation technology.K = 16 users are evenly distributed 7 Comparison of MSE of the proposed angle estimation with the number of Illustrates the variation of angle MSE as SNR for the proposed method with different search times of angle rotation technology in the base station service range, and users are divided into 4 groups according to the user's distribution, and there are 4 users in each group, and L = 16 .The number of searches for the rotation angle is set to N ro = 10 , N ro = 20 , N ro = 40 , N ro = 60 , N ro = 100 and N ro = 200 , respectively.The performance metric of the angle estima- tion is taken as the average individual angle MSE, i.e., It can be seen from Fig. 6 that the performance of angle estimation is better as the number of searches increases.However, when the number of searches increases to a certain value and continues to increase the number of searches, the performance improvement for angle estimation is not very large.Therefore, the number of searches in proposed method is set to N ro = 60.

Conclusions
In this paper, the channel estimation problem of single-cell UAV 3D MIMO systems is studied.As the number of antennas increases, the conventional channel estimation methods (such as 2D-DFT estimation method) will greatly increase the computational complexity of the systems.Moreover, as the number of users or the number of antennas of the user in the cell increases, the pilot overhead of the system also increases.In order to solve the problems of computational complexity and pilot overhead of the system, this paper proposes a low complexity channel estimation algorithm using paired spatial signatures.Firstly, 3D channel modeling based on the angle between the direction of arrival and the x-direction of the array antenna is projected onto the x-and y-directions to perform channel estimation separately respectively, followed by the direction of arrival and the y-direction of the array antenna, and then the 3D-MIMO channel.Finally, the 3D MIMO channel is generated by the Kronecker product.In addition, this paper also proposes a 3D-ADMA group and a two-stage pilot allocation strategy to reduce the pilot overhead of the system.Compared with the conventional channel estimation method of a 3D MIMO system UPA using a low-rank model, the proposed methods greatly reduce the computational complexity without reducing the estimated performance.For subsequent work, we will consider low complexity 3D MIMO channel estimation methods with different channel gains and multipath conditions.
Figure 2 depicts a system model of users' signals reaching UAV-based BS in a cell.

hFig. 2
Fig. 1 UPA signal model.This is UPA signal model.Depicts a UPA signal model.The antenna is evenly distributed on the xoy plane, and the antenna spacing of the x-axis and the y-axis is d.φ k , θ k denote the signal elevation angle and azimuth of the kth user respectively.α k , β k denote the angle between the incident signal and the x-and y-directions of the kth user respectively

Fig. 3
Fig.3The transmission data phase diagram.Shows the system data transmission, which is divided into two stages: uplink preamble and data transmission.The data transmission stage includes the uplink and downlink stages

Fig. 4 1 σ 1 .
Fig. 4 Number of MFLOPs versus number of antennas, L = 16 G = 4 , g = 4 , N ro = 60 N co = 10 , N b = 10 .Shows the comparison of different algorithm computational complexity.Curve of the number of MFLOPs changes with the number of antennas

Fig. 5
Fig. 5 Comparison MSE performances of the proposed method with L = 16, 32, and 64, respectively.Illustrates the MSE performance as a function of SNR for the proposed methods with lengths of pilot sequences L = 16, L = 32, and L = 64, respectively

Table 1
Computational complexity