I Introduction
Terahertz (THz) communication is aimed at playing a vital role in the future sixth generation (6G) wireless communication networks [1]. To support up to terabit per second (Tbps) ultrahigh peak data rate, THzband, whose spectrum ranges from 0.1 to 10 THz, can offer tens of gigahertz (GHz) ultrabroadband bandwidths [2]. Meanwhile, ultramassive multipleinput multipleoutput (UMMIMO)based transceiver, which is equipped with tens of thousands of antennas [3], would be realized in THz communications. The UMMIMO array of THz transceiver can utilize its large enough array aperture and beamforming techniques to effectively combat the severe path loss of THz signals and further extend the communication range [4]. Therefore, THz UMMIMO technique is a favorable candidate for the unmanned aerial vehicles (UAVs) and satellite communications [4].
The reliable channel state information (CSI) acquisition is indispensable for the qualityofservice (QoS) of this spacetoair scenario [5, 6, 7]. However, the highspeed mobility of UAVs and satellites makes accurate CSI acquisition rather challenging. To acquire the dominant accurate channel parameters including the angles, Doppler shifts, and channel gains, some multistage channel estimation solutions were proposed in [8] for narrowband millimeterwave (mmWave) MIMO systems. In [9], a prioriaided THz channel estimation scheme was proposed to predict and track the physical direction of lineofsight (LoS) component of the timevarying massive MIMO channels in THz beamspace domain.
However, due to the unprecedentedly ultrabroad band, ultralarge array aperture, and ultrahigh velocity in the THz UMMIMO based satellite communication systems, the aforementioned channel estimation schemes are difficult to be applied. Compared with the mmWave MIMO systems with limited bandwidth and aperture, the THz UMMIMO based spacetoair channels present the unique triple delaybeamDoppler squint effects. Specifically, the THz UMMIMO array equipped with tens of hundreds of antennas will suffer from different propagation delays at different antennas for the same impinging signal. These delay gaps may be as large as several symbol periods in ultrabroadband THz communications, which indicates the nonnegligible intersymbolinterference even for the LoS link. We can call this inevitable phenomenon in THz UMMIMO systems as the delay squint effect, which is also named as spatialfrequency wideband effects in [10, 11]. Meanwhile, the beam squint effect, in which the beam direction is a function of the operating frequency, can be further introduced by this delay squint effect. The large Doppler shift caused by the highspeed mobility in spacetoair THz communications is also frequencydependent, i.e., Doppler squint effect. The THz UMMIMO based UAVs and satellite communication systems present triple delaybeamDoppler squint effects. Nevertheless, recent researches mainly focus on the impact of beam squint effect on mmWave systems [10, 11, 12]. Consequently, an efficient channel parameter estimation solution is indispensable for THz UMMIMO based spacetoair communications.
In this paper, we mainly investigate the angle estimation for spacetoair LoS links connecting multiple UAVs and low earth orbit (LEO) satellite^{1}^{1}1Due to space limitation, we only estimate the azimuth and elevation angles at UAVs and satellite, while the remaining parts of channel estimation and tracking can be found in our work published in IEEE JSAC [13].. To save space, we assume Doppler shifts can be well compensated using the positioning and flight posture information acquired by satellite communications, so that we just consider the dual delaybeam squint effects of spacetoair THz UMMIMO channel (without Doppler squint effect) in this paper. Specifically, based on the rough angle estimates acquired from navigation information and positioning system, we first design a grouping truetime delay unit (GTTDU) module with low hardware cost to significantly mitigate the impact of delaybeam squint effects on both the transmitter and receiver, so as to establish the spacetoair THz link. After the link establishment, the UM hybrid array can be equivalently considered as a lowdimensional fullydigital array based on the subarray selection scheme. Then, the fine estimates of azimuth/elevation angles at both UAVs and satellite can be separately acquired using the proposed prioraided iterative angle estimation algorithm. Simulations results have the good tightness with the CramérRao lower bounds (CRLBs) of azimuth/elevation angles, which testifies the good performance of the proposed solution.
Notations
: Italic boldface lower and uppercase symbols denote column vectors and matrices, respectively.
, , and denote the conjugate, transpose, and Hermitian transpose operations, respectively. is the norm of . and denote the Kronecker and Hadamard product operations, respectively. denotes the vector of size with all the elements being . is the cardinality of the set , and denotes the th element of the ordered set . denotes the subvector containing the elements of indexed in the ordered set . and denotes the th element of and the throw and the thcolumn element of , respectively. Finally, is the expectation of the argument.Ii THz UMMIMO Channel Model
In this section, we will formulate the THz fulldimensional UMMIMO channel model using uniform planar array (UPA), which involves azimuth and elevation angles [14]. Fig. 1(a) depicts the specific scenario that UAVs communicate with a LEO satellite through respective THz LoS link. In Fig. 1(b), we can observe that the hybrid beamforming structure with a subconnected phase shift network (PSN) can be adopted in the transceiver at LEO satellite and UAVs [2]. The specific configurations of these antenna arrays are as follows. Defining and as the numbers of antennas in horizontal and vertical directions, respectively, and the total number of antennas at UAV arrays is . Considering the subconnected PSN at LEO satellite, () and () can be defined as the numbers of subarrays (antennas within each subarray) in horizontal and vertical directions, respectively; while the numbers of antennas in horizontal and vertical directions of array are and , respectively. Then, the total numbers of antennas in the whole antenna array (each subarray) is (). Clearly, RF chains and only one RF chain are equipped at the LEO satellite and UAV, respectively, and each RF chain and the corresponding subarray mounted on LEO satellite are allocated to the specific UAV.
To combat the multipath effect at the receiver of LEO satellite caused by multiple THz LoS links, the orthogonal frequencydivision multiplexing (OFDM) technique with subcarriers will be applied to this spacetoair communication system. Specifically, the orthogonal frequency division multiple access (OFDMA) is utilized to split the pilot signals transmitted from different UAVs and improve the accuracy of angle estimation. Hence, considering the alternating subcarrier index allocation with equal intervals, subcarriers can be equally assigned to UAVs, and the ordered subcarrier index set corresponding to the th UAV is with . Due to the channel reciprocity of downlink (DL) and uplink (UL) in time division duplex (TDD) systems, we focus on the formulation of DL channel matrix next. The DL spatialdelay channel matrix at time corresponding to the th UAV can be defined as , whose the ()th element, i.e., , can be expressed as
where , , is the largescale fading gain, denotes the channel gain^{2}^{2}2Due to the negligible attenuation of THz spacetoair communication links (e.g., atmospheric molecular absorption) in the stratosphere and above, the channel gain can be modeled as a frequency flat coefficient., is the Doppler shift with and being the carrier wavelength at central carrier frequency and the relative radial velocity, respectively, is the path delay, and denote the transmission delay for the th and th antennas at the UAV and satellite arrays, respectively, and is the Dirac impulse function. After some algebraic transformations, the DL channel matrix
in spatialfrequency domain at the
th subcarrier of the th OFDM symbol for the th UAV can be formulated as(1) 
where and are the system bandwidth and the duration time of an OFDM symbol, respectively, and denotes the DL array response matrix associated with the th UAV, given by
(2) 
where () and () that consider halfwavelength antenna spacing are the virtual angles at horizontal and vertical directions for the th UAV (satellite), respectively, is the general DL array response matrix without beam squint effect, and is the array response squint matrix with beam squint effect. In (2), and are the regular array response vectors at the th UAV and satellite [13, 14], respectively, and and are the frequencydependent array response squint vectors, respectively. Moreover, the expressions of the horizontal/vertical steering vectors and , and the horizontal/vertical steering squint vectors and can be found in [13] for details. Note that the vectors at UAVs, i.e., , , , and , have the similar definitions and expressions.
Similar to (II), the UL spatialfrequency channel matrix at the th subcarrier of the th OFDM symbol for the satellite, denoted by , can be expressed as
(3) 
Iii Proposed Angle Estimation Solution under DelayBeam Squint Effects
In this section, we propose a prioraided iterative angle estimation solution for THz UMMIMO based spacetoair communications, where some rough prior angles at UAVs and LEO satellite can be acquired from the positioning, flight speed and direction, and posture information for establishing the THz communication links.
To overcome the delaybeam squint effects of THz UMMIMO array, we design a hardwareeffective implementation of TTDU module, i.e., GTTDU module based transceiver structure as shown in Fig. 2, as the alternative of optimal TTDU module that is made up of numerous truetime delay units with each unit being assigned to its dedicated antenna [15]. From Fig. 2, we observe that except for the antenna array, this transceiver structure contains a GTTDU module and a reconfigurable RF selection network involving a subconnected PSN and an antenna switching network (ASN), where this ASN can control the active or inactive state of the antenna elements to form different connection patterns of the RF selection network at the angle estimation stage. In this GTTDU module, a TTDU can be shared by a group of antennas and this imperfect hardware limitation can be handled by the subsequent signal processing algorithms well.
Based on the prior information acquired from positioning systems, the rough estimates of azimuth and elevation angles at UAVs (satellite) can be defined as () and (), respectively, and the corresponding horizontally and vertically virtual angles are () and (), respectively. According to in (II), we present the DL spatialfrequency channel matrix compensated by ideal TTDU module as
(4) 
in which
(5) 
By comparing in (5) and in (2), we can find that the beam squint effect part can be perfectly eliminated if we can acquire the perfect angle information, i.e., and then when , , , and . Moreover, the compensated UL spatialfrequency channel matrix has the similar expressions, which are omitted for simplicity. The practical DL/UL channel matrices compensated by the GTTDU module can be derived from (III) and (5).
At the angle estimation stage, the azimuth/elevation angles at UAVs can be estimated in DL, while those angles at satellite are acquired in UL.
Iiia Fine Angle Estimation at UAVs
Due to the limited valid observation at the UAVs, it is necessary to accumulate multiple OFDM symbols for estimating the azimuth and elevation angles. Meanwhile, the transmitted signals can be compensated the priori Doppler shifts well, that is, the compensated channels will be slow timevarying. By altering the RF connection pattern at the antenna array of UAV, the received signals that adopt different selected subarrays will differ by one envisaged phase value if the transceiver has the same configuration, and these organized phase differences can construct the array response vector of lowdimensional fullydigital array. In Fig. 3, we take the UPA of size as an example. By controlling the reconfigurable RF selection network subarrays of size in successive OFDM symbols can be selected to form the array response vector of equivalent fullydigital array with size of with the critical antenna spacing .
Specifically, OFDM symbols are used to estimate the angles at UAVs, where each OFDM symbol adopts a selected subarray (corresponding to the dedicated RF connection pattern). By employing the rough angle estimates at satellite and UAVs, we can first design the analog precoding and combining vectors, i.e., and for , . For , initialize as , and then let , where with denotes the antenna index of subarray assigned to the th UAV due to each subarray at satellite only communicating with its corresponding UAV as shown in Fig. 1(b). To design , we can partition the UMMIMO array at UAV into smaller subarrays for yielding the array response vector of equivalent lowdimensional fullydigital array with size of . The sizes of these smaller subarrays are ( and ) and their number of antennas is . For and , by defining with and being the ()th subarray, respectively, we denote the antenna index of the selected th subarray corresponding to the th OFDM symbol as , where So can be also initialized as , and then let for .
The received signal at the th subcarrier of the th OFDM symbol transmitted by the th UAV is given by
(6) 
where , , is the channel matrix compensated by the Doppler shifts and the GTTDU module, and and are the transmitted pilot signal and noise, respectively. By collecting the received signals at subcarriers as , we have
(7) 
where , is the error vector including the residual beam squint caused by inaccurate prior information. Moreover, the same transmitted pilot signals are adopted for OFDM symbol, i.e., for . By stacking received from OFDM symbols, we can obtain , i.e.,
(8) 
where and are the residual beam squint and the analog combining matrices, respectively. Based on this analog combining matrix , we can extract the array response vector of equivalent lowdimensional fullydigital to estimate the angles at UAVs using the robust array signal processing techniques. Specifically, compared with for in (7), is multiplied by an extra phase shift for and . Obviously, the effective array response vector of equivalent fullydigital array with size of can be constituted by these regular phase shifts, that is, . Thus, in (IIIA) can be then rewritten as
(9) 
where denotes the beamaligned effective channel gain.
Next, a prioraided iterative angle estimation is proposed as follows. By applying the twodimensional unitary ESPRIT (TDUESPRIT) algorithm [14] to the received signal matrix in (9), the proposed algorithm at the first iteration estimates the azimuth and elevation angles (the corresponding horizontally and vertically virtual angles) at the th UAV, i.e., and ( and ) for . Furthermore, to minimize the impact of on (9), the estimated angles above can be utilized to iteratively compensate at the subsequent iterations (i.e., ), so as to acquire more accurate angle estimates. To be specific, according to the roughly priori virtual angle estimates and , and and estimated at the th iteration, the compensation matrix at the th iteration can be defined as , whose the th column can be expressed as
(10) 
The processed matrix compensated by is
(11) 
To obtain the more accurate angle estimates, we can apply the TDUESPRIT algorithm to those matrices until the maximum iterations is reached, i.e., and the estimated azimuth/elevation angles at UAVs and the corresponding virtual angles are , , , and for .
Algorithm 1 summarizes the proposed prioraided iterative angle estimation procedure above, which can address the beam squint effect well.
IiiB Fine Angle Estimation at Satellite
According to the TDD channel reciprocity, estimating the fine angle estimates at satellite in UL is similar to the DL angle estimation of UAVs in Section IIIA
. At this stage, the fine angles estimated at UAVs above can replace the originally rough angle estimates to refine the GTTDU modules at UAVs and design the analog precoding vectors at UAVs for beam alignment with improved receive signaltonoise ratio (SNR).
To be specific, considering OFDM symbols to estimate the fine azimuth/elevation angles at satellite, we can obtain the lowdimensional fullydigital array with size of . By employing the estimated , we design the analog precoding vector as for . According to the reconfigurable RF selection network, we denote the selected antenna index in the th OFDM sysmbol at the th aircraft subarray as , where . The analog combining vector can be then initialized as , and let , for , .
Similar to (6)(9), the received signal can be also obtained, where the details can be found in [13] due to the limited space. By replacing the input parameters for UAVs with the corresponding parameters for satellite, we can utilize the proposed prioraided iterative angle estimation algorithm in Algorithm 1 to acquire more accurate estimates of azimuth and elevation angles and the corresponding virtual angles at satellite, i.e., , , , and for .
IiiC Computational Complexity Analysis
The computational complexity of the proposed angle estimation solution consists of the acquisition of azimuth/elevation angles at UAVs and satellite using TDUESPRIT algorithm. Their total computational complexity is . Since the effective lowdimensional signals at the UAVs and satellite are utilized to estimate the angles, the computational complexity of the proposed angle estimation solution is in polynomial time even although the THz UMMIMO arrays are equipped with tens of thousands of antennas at UAVs and satellite.
Iv Numerical Evaluation
In this section, we evaluate the performance of the proposed angle estimation solution for THz UMMIMObased spacetoair communications. Without loss of generality, the LEO satellite serves UAVs, where the vertical distance between satellite and UAVs is kilometer and these two UAVs are randomly appeared in a horizontal circular plane with radius . The relative radial velocity among UAVs and satellite is meter per second. Other simulation parameter settings are shown in Table I. Moreover, we define with
being the noise variance as the transmitted SNR in DL and UL. The performance of angle estimation is evaluated using the root mean square error (RMSE) metric given by
, where and represent the true and the estimated angle vectors with being , , , or . To evaluate the estimation performance, the CRLBs serve as the lower bounds of angle estimation [16].


Parameter  Value 


()  THz ( GHz) 
, , , 