1 Introduction↩︎

With the successful deployment of the fifth-generation (5G) mobile communication systems worldwide, research on the sixth-generation (6G) wireless networks has intensified[1]. To support the ambitious objectives of 6G wireless networks[[2]][3], such as unprecedented high data rates, ultra-high reliability, global coverage, and ultra-dense connectivity, a variety of new solutions and technologies have recently emerged. Among these, integrated sensing and communications (ISAC)[4], reconfigurable intelligent surfaces (RIS)[5], and extremely large antenna array (ELAA)[[6]][7], also known as extremely large-scale array (XL-array) [8] are being actively investigated as potential enabling technologies for 6G systems.

Despite their potential to enhance communication performance, these technologies face challenges such as severe free-space path loss and sensitivity to blockage. To achieve ubiquitous 6G coverage, unmanned aerial vehicles (UAVs) have been widely studied for their mobility, flexibility, and adaptive altitude adjustment [9] [10]. When deployed as aerial base stations (BSs), UAVs enhance wireless coverage and establish favorable propagation conditions, such as line-of-sight (LoS) distribution. However, due to limited payload capacity and LoS-dominated air-to-ground links, a single UAV-enabled antenna array cannot ensure scalable aperture gain and spatial multiplexing [11]. Multi-UAV systems address these challenges by extending coverage, ensuring reliable LoS transmission, and improving task efficiency through collaboration while maintaining maneuverability [12]

Among the numerous envisioned 6G enabling technologies, ELAA has emerged as a promising solution to meet the ultra-high spectral efficiency and spatial resolution demands of future networks [[6]][13]–[17]. However, the significant increase in antenna size poses deployment challenges, especially for mobile platforms like vehicles, aircraft, and ships, where space for array antennas is fragmented and limited. To address this, researchers have explored distributed ELAA, including modular ELAA [18] and widely-spaced multi-subarray (WSMS) concepts [[19]][20], which enlarge the antenna aperture without additional elements. By dividing a large-aperture array into smaller subarrays deployed across mobile platforms, ELAA performance is maintained while adhering to mobility constraints.

Given this context, we propose to investigate the multi-UAV near-field system, as the increased array aperture brought by multi-UAV makes it easier to enter the near-field, which will present several advantages. First, multi-UAV near-field systems can expand the overall antenna aperture without additional elements, thereby achieving high spatial resolution and enabling precise positioning. Secondly, by equipping each UAV with a medium-sized array, distributed configurations enhance both flight stability and operational flexibility. Finally, this integration complements non-terrestrial networks, fostering the advancement of future wireless communications and their applications.

1.1 Related Work↩︎

1.1.1 Near-field performance analyses↩︎

The deployment of ELAA fundamentally alters channel characteristics, transitioning from far-field to near-field communication, where the conventional plane wave model (PWM) becomes inaccurate. Instead, the spherical wave model (SWM) better captures phase and amplitude variations across large arrays[14], [15], [18]–[20]. Prior works have derived closed-form SNR expressions and asymptotic performance analyses under SWM, revealing that SNR scaling deviates significantly from the linear trend predicted by PWM, highlighting the necessity of near-field modeling.

For the distributed ELAA, it is not only capable of accommodating a very large number of antennas, but also meets the practical constraints of the installation structure. More importantly, as the spacing increases, it becomes easier to access the near-field region, which allows for the full exploitation of the advantages inherent to the near-field environment. To accommodate its spatial structure, the hybrid spherical and plane wave model (HSPWM) is adopted, where the planar wavefront is applied to the intra-subarray and the spherical wavefront is used for the inter-subarray[[19]][20][[21]][22]. HSPWM is also referred to as a cross-field model, which captures the transition between near- and far-field regimes. [19] derived the closed-form angle and range Cramér-Rao bounds (CRBs) for the spherical wavefront-based WSMS (SW-WSMS) and HSPWM-based WSMS (HSPWM-WSMS) in bi-static systems. The potential for high-resolution near-field localization using WSMS has been verified through comparative analysis. On the other hand, the application of the HSPWM assumption not only circumvents the inaccuracies associated with the PWM but also mitigates the increased parameterization that arises from the expanded channel matrix dimensions during channel estimation when employing the SWM [[21]][22]. Additionally, it also enriches the multiplexing within the propagation paths[20]. Additionally, by analyzing both near-field and far-field channel models, [23] proposes a two-stage hybrid-field beam training scheme based on coarse beam sweeping and fine hierarchical codebooks. This approach distinguishes near-field users (NUs) from far-field users (FUs) using a modified Rayleigh distance (defined by 90% beam gain) and 3dB spatial angle width.

1.1.2 Near field channel estimation↩︎

Accurate Channel State Information (CSI) is vital for ELAA, but conventional far-field estimation methods degrade significantly due to the near-field effect. Recent works have proposed dedicated near-field solutions [24]–[26]. These include methods based on joint angle-distance sparse representations [25] and low-complexity approaches that decouple these parameters for UPAs [26].

All of these methods are based on uniform ELAA, where antenna elements are deployed with uniform inter-element spacing, such as in uniform linear or planar arrays. In contrast, research has also been conducted on channel estimation methods for distributed ELAA, where subarrays are geographically separated. Literature [21] divided array antennas into smaller subarrays, treating each subarray’s channel as a far-field channel. A two-stage channel estimation scheme was used, with a deep convolutional neural network estimating subarray channel parameters in the first stage and geometric relationships completing the estimation in the second. To simplify complexity, the HSPWM assumption was employed [[19]][27], enabling a two-dimensional orthogonal matching pursuit (2D-OMP) method with lower complexity than traditional polar domain orthogonal matching pursuit (PD-OMP) [27], reducing the computational burden for practical applications.

Near-field channel estimation becomes particularly critical in multi-UAV-assisted communication systems due to the limited altitude, rapid topology changes, and short communication distances, which frequently result in near-field propagation conditions. Several recent studies have applied deep reinforcement learning and generative models to optimize multi-UAV beamforming and secure communications in complex environments [28], [29]. These works emphasize the growing importance of robust and adaptive channel estimation techniques under near-field and distributed multi-UAV scenarios.

1.2 Contributions↩︎

Although existing literature has extensively explored subarray-based channel modeling, the integration of distributed characteristics and near-field propagation conditions in multi-UAV communication systems remains underexplored. Given that individual UAV is typically equipped with UPA, it is of significant importance to investigate the performance of multi-UAV near-field communication systems. Moreover, existing UAV channel estimation methods are mostly based on far-field assumptions, primarily due to the limited number of antennas on a single UAV and the considerable distance from the ground, which rarely satisfy near-field conditions. However, due to the increased aperture caused by UAV deployment scale and spacing, multi-UAV scenarios are likely to operate in the near-field region. For example, in a \(2 \times 2\) UAV configuration with inter-UAV spacing of \(50d\) (\(d\) denotes the half-wavelength), the Rayleigh distance exceeds 110 m at 10 GHz. This motivates our research on near-field channel estimation for multi-UAV systems. Consequently, this work establishes the first comprehensive framework for analyzing both performance metrics and channel estimation in multi-UAV near-field systems. To our knowledge, it pioneers the application of ELAA principles to practical UAV swarm configurations. We have made the following contributions:

• We model a multi-UAV to user communication framework in which each UAV is equipped with a UPA. The UAV spacing is larger than half-wavelength to cooperatively perceive the user in the ground. To ensure general applicability, we adopt a rectangular formation for the UAV swarm configuration.

• We investigate the mathematical modeling and SNR performance analysis in multi-UAV near-field communication system. The analytical expressions for the maximum SNR under three models including SWM, PWM and HSPWM are first derived. Performance analysis reveals the unique potential of HSPWM in distributed ELAA systems, characterized by its simplicity and high precision. Notably, this is the first work to conduct performance analysis under the assumption of HSPWM for the distributed ELAA and to compare it comprehensively with both PWM and SWM.

• We propose a channel estimation method applying the HSPWM assumption. First, we extend the PD-OMP method directly to the multi-UAV near-field scenario, thereby developing the spherical domain orthogonal matching pursuit (SD-OMP) method. However, it does not fully exploit the structural characteristics of the near-field array response, leading to a large dictionary size and consequently high complexity in channel estimation. To address this, we further propose a tensor dictionary based on HSPWM, facilitating the development of a low-complexity near-field channel estimation method. Moreover, we employ the Nelder-Mead algorithm to resolve the mismatch issue.

The paper is organized as follows. Section II introduces the system model of the multi-UAV distributed near-field system. Section III analyzes the SNR performance. Channel estimation methods are proposed in Section IV. Simulation results are shown in Section V. Section VI concludes the paper.

Notations: \(a\), \(\boldsymbol{a}\), \(\boldsymbol{A}\) and \(\boldsymbol{\mathcal{A}}\) represent the scalar, the vector, the matrix and the tensor, respectively. \(\mathbb{C}^{M \times N}\) depicts the set of \(M \times N\)-dimensional complex-valued matrices. \((\cdot)^{T}\), \((\cdot)^{H}\) and \((\cdot)^{{-1}}\) represent transpose, conjugate transpose and inverse operation, respectively. \((\cdot)^{{\dagger}}\) denotes the pseudoinverse, and its formula is \(\boldsymbol{A}^{\dagger} = (\boldsymbol{A}^{\mathrm{H}} \boldsymbol{A})^{-1} \boldsymbol{A}^{\mathrm{H}}\). \(\exp\{\cdot\}\) defines the exponential function of \(e\). \(\mathbf{I}_N\) defines an \(N\)-dimensional identity matrix. \(\sum_{n=0}^{N}(\cdot)\) represents the sum of all function values from 0 to \(N\). \(\circ\) represents vector outer product. \(\otimes\) represents Kronecker product. \(\times_1\), \(\times_2\) and \(\times_3\) represent the 1-mode, 2-mode and 3-mode product of a tensor with a vector, respectively. \(\|\cdot\|\) and \(\|\cdot\|_F\) stands for the 2-norm and the Frobenius norm. \(\mathcal{CN}(\boldsymbol{a}, \boldsymbol{A})\) denotes the Gaussian distribution in the complex domain, where \(\boldsymbol{a}\) represents the mean, and \(\boldsymbol{A}\) signifies the covariance matrix. \(\mathcal{U}\left(-a, a\right)\) signifies a uniform distribution across the interval \((-a,a)\).

2 SYSTEM MODEL↩︎

As shown in Fig. 1, we consider a near-field SIMO wireless communication system based on uplink time division duplexing (TDD), which comprises \(M=M_x\times M_y\) UAVs and a single-antenna user. Each UAV is equipped with a UPA of \(N=N_x\times N_y\) antennas mounted at its bottom, which is shown in Fig. 1. Specifically, the UAV located in the \(m_x\)-th row and \(m_y\)-th column is indexed by \((m_x, m_y)\), \(m_x = 1, 2, \cdots, M_{x}\) and \(m_y = 1, 2, \cdots, M_{y}\). The index of the antenna in the \(n_x\)-th row and \(n_y\)-th column of a UAV is denoted as \((n_x, n_y)\), \(n_x = 1, 2, \cdots, N_x\) and \(n_y = 1, 2, \cdots, N_y\). For simplicity in notation, in the following paper, the UAV positioned at the \(m_x\)-th row and \(m_y\)-th column is referred to as \(m\), where \(m = (m_x-1) M_x + m_y\); each antenna element on the UAV at the \(n_x\)-th row and \(n_y\)-th column is denoted by \(n\), where \(n = (n_x-1) N_x + n_y\), thereby ensuring the uniqueness of the mapping relationship. Each UAV carries a radio frequency (RF) chain, with each linked to \(\textit{N}\) antenna elements. Therefore, the total number of array elements is \(\textit{MN}\). For analytical simplicity, we assume that all UAVs are perfectly synchronized in time and frequency and remain in a stable hovering state with no relative motion. ²

All communication nodes are placed in the three-dimensional (3D) Cartesian coordinate system. The UAV located in the 1-th row and 1-th column, denoted by \(\textit{U}_1\), serves as the reference. Note that UAVs should not ascend and descend frequently to avoid energy consumption, so we assume that the UAV swarm flies at a specific altitude of \(\textit{H}\).

For each UAV, the first antenna is chosen as the reference antenna. The antenna spacing \(\textit{d}=\lambda/2\) within a UAV, where \(\lambda\) denotes the carrier wavelength. The UAV spacing is arbitrarily multiple of half-wavelength, which represents the general setting of the distributed ELAA. Therefore, the distances of the reference antenna element between adjacent UAVs alone on the \(x\) axis and the \(y\) axis are represented as \(D_{x}d\) and \(D_{y}d\), respectively, where \(D_{x}\), \(D_{y} \gg 1\). So the relative position of the \((m_x, m_y)\)-th UAV with reference to \(\textit{U}_1\) is denoted by \(\boldsymbol{d}_{m}=[(m_{x}-1)D_{x}d, (m_{y}-1)D_{y}d, 0]^T\). Denoting the location of \(\textit{U}_1\) as \(\boldsymbol{p}_{U_1}\), the location of the \((n_x, n_y)\)-th antenna on the \((m_x, m_y)\)-th UAV is \(\boldsymbol{p}_{m, n}=\boldsymbol{p}_{U_1}+\boldsymbol{d}_{m}+\boldsymbol{d}_{n}\), where \(\boldsymbol{d}_{n}=\left[(n_{x}-1)d, (n_{y}-1)d, 0\right]^T\) is the position of the antenna relative to the first antenna of that UAV.

The user’s location is denoted as \(\boldsymbol{p}_E\), and its relative position to \(\textit{U}_1\) is \(\boldsymbol{p}_{U_1E} = \boldsymbol{p}_{U_1} - \boldsymbol{p}_{E} =\left[ r \Psi, r \Phi, r \Omega \right]^T\), where \(r=\left\| \boldsymbol{p}_{U_1} - \boldsymbol{p}_{E} \right\|\), \(\Psi=\sin \varphi \cos \theta\), \(\Phi=\sin \varphi \sin \theta\), and \(\Omega=\cos \varphi\), with \(\theta \in [0, \pi]\) and \(\varphi \in \left[-\frac{\pi}{2}, \frac{\pi}{2}\right]\) denote the azimuth and the elevation angle, respectively. Thus the distance between the user and the \((n_x,n_y)\)-th antenna element of the \((m_x,m_y)\)-th UAV is then given by

where \(\xi = \frac{d}{r}\ll 1\), and \(r=r_{1,1}\) is the distance between the user and the reference antenna of \(\textit{U}_1\).

2.1 Signal model↩︎

Assuming that the pilot signal sent by the user at \(q\)-th time slot is denoted as \(x_q\), the received signal \(y_{q,m}\) on the \(m\)-th UAV can be expressed as:

where \(\boldsymbol{w}_{q,m}\in \mathbb{C}^{N \times 1}\) represents the measurement combiner at the UAV side , and \(\boldsymbol{h}_{m} \in \mathbb{C}^{N\times1}\) denotes the wireless channel between the \(m\)-UAV and the user. \(\boldsymbol{n}_{q,m}\) denotes the complex additive white Gaussian noise (AWGN), which follows the distribution \(\mathcal{CN}(0,\sigma^2\mathbf{I}_N)\).

Assuming that the pilot signal \({x}_q\) sent by the user is 1 and \(Q\) is denoted as the pilot length, the collection of all received signals \(\boldsymbol{y}_m\) on the \(m\)-th UAV can be expressed as:

where \(\boldsymbol{W}_m = [\boldsymbol{w}_{1,m}, \ldots, \boldsymbol{w}_{Q,m}] \in \mathbb{C}^{N \times Q}\) denotes the measurement matrix of the \(m\)-th UAV subarray. For a UPA with \(N=N_x N_y\), the measurement matrix \(\boldsymbol{W}_m\) is designed as \(\boldsymbol{W}_m=\boldsymbol{W}_x \otimes \boldsymbol{W}_y\), with \(\boldsymbol{W}_x \in\mathbb{C}^{N_x \times Q_x}\) and \(\boldsymbol{W}_y \in\mathbb{C}^{N_y \times Q_y}\) representing the measurement matrix along the \(x\)-axis and \(y\)-axis, respectively. \(Q_x\) and \(Q_y\) denote the number of measurements along the \(x\)-axis and \(y\)-axis, respectively. Additionally, \(\overline{\boldsymbol{n}}_m = [\boldsymbol{w}_{1,m}^H \boldsymbol{n}_{1,m}, \ldots, \boldsymbol{w}_{Q,m}^H \boldsymbol{n}_{Q,m}]^T \in \mathbb{C}^{Q \times 1}\) represents the noise projection.

According to compressed sensing theory, to enable efficient channel estimation, the measurement matrices \(\boldsymbol{W}_x\) and \(\boldsymbol{W}_y\) are required to exhibit low mutual coherence. Taking \(\boldsymbol{W}_x\) as an example (the construction of \(\boldsymbol{W}_y\) follows similarly), we construct \(\boldsymbol{W}_x \in \mathbb{C}^{N_x \times Q_x}\) based on singular value decomposition (SVD) as:

where \(\boldsymbol{U}_x \in \mathbb{C}^{N_x \times N_x}\) and \(\boldsymbol{V}_x \in \mathbb{C}^{Q_x \times Q_x}\) are unitary matrices, and the singular value matrix is given by \(\boldsymbol{\Sigma}_x = \begin{bmatrix} \mathrm{diag}(\boldsymbol{\sigma}) , \mathbf{0}_{Q_x,\, N_x - Q_x} \end{bmatrix}^T\), and\(\quad \boldsymbol{\sigma} = [\sigma_1, \ldots, \sigma_{Q_x}]^T\) and \(\mathbf{0}_{Q_x,\, N_x - Q_x}\in \mathbb{C}^{Q_x \times (N_x-Q_x)}\) is a null matrix.

To ensure orthogonality, the following Frobenius norm metric is introduced:

This expression is minimized when all singular values equal one, i.e., \(\sigma_1 = \cdots = \sigma_{Q_x} = 1\). Accordingly, the unconstrained optimal solution is given by:

where \(\boldsymbol{U}_x \in \mathbb{C}^{N_x \times N_x}\) and \(\boldsymbol{V}_x \in \mathbb{C}^{Q_x \times Q_x}\) are arbitrary unitary matrices.

To meet the modulus-1 constraint imposed by phased-array hardware, each element of \(\boldsymbol{W}_x\) is normalized as:

By collecting signals from all UAVs, the received signal \(\boldsymbol{y}\) can be represented as:

where \(\boldsymbol{y} = \begin{bmatrix} \boldsymbol{y}_1^T, \ldots, \boldsymbol{y}_M^T \end{bmatrix}^T \in \mathbb{C}^{MQ \times 1}\), \(\boldsymbol{h} = \begin{bmatrix} \boldsymbol{h}_1^T, \ldots, \boldsymbol{h}_M^T \end{bmatrix}^T \in \mathbb{C}^{MN \times 1}\) and \(\overline{\boldsymbol{n}} = \begin{bmatrix} \overline{\boldsymbol{n}}_1^T, \ldots, \overline{\boldsymbol{n}}_M^T \end{bmatrix} \in \mathbb{C}^{MQ \times 1}\). Here, \(\boldsymbol{W}\in \mathbb{C}^{MN \times MQ}\) denotes the overall measurement matrix composed of \(M\) subarrays. In this study, we assume that each subarray employs the same measurement matrix, i.e., \(\boldsymbol{W}_0=\boldsymbol{W}_1=\ldots=\boldsymbol{W}_{M-1}\). Consequently, \(\boldsymbol{W}\) has a block diagonal structure and can be streamlined to \(\boldsymbol{W}=\mathbf{I}_{M} \otimes \boldsymbol{W}_x \otimes \boldsymbol{W}_y\), which can be expressed as

2.2 Channel model↩︎

Considering the geometric structure and the nature of limited scattering, the channel \(\boldsymbol{h}\) is written as

where \(\textit{L}\) is the number of channel paths, \(\textit{z}_l\) denote the complex gain of the \(\textit{l}\)-th path. As depicted in Fig. 1, \(\{\theta_l\}_{l=1}^L\), \(\{\varphi_l\}_{l=1}^L\) and \(\{r_l\}_{l=1}^L\) denote the azimuth angle, the elevation angle and distance parameters, respectively. The array response \(\boldsymbol{g}(\theta, \varphi, r)\) is given by \(\boldsymbol{g}(\theta, \varphi, r)=[g_{0,0}(\theta, \varphi, r),\ldots,g_{m,n}(\theta, \varphi, r),\ldots,g_{M,N}(\theta, \varphi, r)]^T\in \mathbb{C}^{MN \times 1}\), which \(g_{m, n}\left(\theta, \varphi, r\right)\) is the array response element between the user and the \((n_x, n_y)\)-th element antenna in the \((m_x, m_y)\)-th UAV. \(g_{m, n}\left(\theta, \varphi, r\right)\) can be expressed as \({g}_{m, n}(\theta, \varphi, r) = \frac{\sqrt{\beta_{0}}}{r_{m, n}} e^{-j \frac{2\pi}{\lambda} r_{m, n}}\) with \(\beta_{0}\) denoting the channel power at the reference distance \(\textit{d}_0\) = 1 m.

In this subsection, various forms of channel responses tailored to the unique signal propagation characteristics of different communication scenarios are presented, namely \(\boldsymbol{g}_{\mathbf{SWM}}(\theta, \varphi, r)\), \(\boldsymbol{g}_{\mathbf{PWM}}(\theta, \varphi, r)\) and \(\boldsymbol{g}_{\mathbf{HPSW}}(\theta, \varphi, r)\).

2.2.1 SWM↩︎

For users in the near-field region, the channel modeling should follow the SWM, which is the most accurate model that calculates the channel response from all antennas to the user. Consequently, the array response vector \(\boldsymbol{g}_{\mathbf{SWM}}\left(\theta, \varphi, r\right) \in \mathbb{C}^{M N \times 1}\) between the user and the distributed ELAA can be given by (2), in which, the term \(\frac{r}{r_{m, n}}\) determines the magnitude, and the term \(r_{m,n}-r\) determines the phase.

None

2.2.2 PWM↩︎

For far-field scenarios where the user-to-array distance surpasses the Rayleigh distance, the channel model is well approximated by the PWM, given the wavefront’s near-planar nature. Therefore, the amplitude variation caused by \(r_{m,n}\) is assumed to be relatively small. In line with most papers investigating far-field scenarios[[9]][32][33], \(\frac{r}{r_{m, n}}\) is set to 1 in this paper. Additionally, since \(r_{m,n}-r\) determines the phase, even a subtle change will have great impact, which can be represented as

Consequently, the array response vector \(\boldsymbol{g}_{\mathbf{PWM}}\left(\theta, \varphi, r\right) \in \mathbb{C}^{M N \times 1}\) for PWM can be expressed as (3).

It is important to recognize that, unlike the SWM, which depends on both amplitude and phase variations, the PWM is solely governed by phase variation. As a result, the PWM serves as an approximation of the SWM when the array size is much smaller than the communication distance. In other words, when the distance between the user and the array is greater than the Rayleigh distance, the SWM degenerates into the PWM.

2.2.3 HSPWM↩︎

Drawing from above analyses, while the SWM provides the most accurate depiction of channel characteristics, this detailed modeling introduces significant complexities in channel estimation, posing many challenges. On the other hand, the PWM, despite its analytical simplicity, incurs a loss in precision, which diminishes its applicability in near-field conditions. Recognizing these limitations, researchers have turned their focus to the HSPWM assumption, which offers a more balanced approach [19]. This assumption is gaining traction for its potential to mitigate the challenges associated with both the SWM and PWM, providing a more viable framework for channel estimation.

None

Specifically, under the HSPWM assumption, the PWM is utilized within each subarray, maintaining accuracy due to the small array size. The SWM is then applied to enhance the modeling precision within these subarrays. In our scenario, PWM is applied to the UPA within one UAV, while SWM is utilized among UAVs. Considering the HSPWM assumption, \(\boldsymbol{g}_{\mathbf{HSPWM}}(\theta, \varphi, r)\in \mathbb{C}^{M N \times 1}\) is obtained by approximating \(\boldsymbol{g}_{\mathbf{SWM}}(\theta, \varphi, r)\)

where \(\boldsymbol{a}(\theta, \varphi) \in \mathbb{C}^{N \times 1}\) denotes the plane-wave array manifold for each UAV, and \(\boldsymbol{b}(\theta,\varphi, r) \in \mathbb{C}^{M \times 1}\) denotes the spherical-wave array manifold for the inter-UAV. Their expressions are as follows:

3 performance analysis of SNR↩︎

In this section, we initially provide a closed-form expression solution for the SNRs of three models, i.e., SWM, PWM and HSPWM. Furthermore, an asymptotic performance analysis is provided to describe the impact of the distribution of UAV array on SNR in our multi-UAV mmwave communication system. Upon comparison, it is evident that the HSPWM delivers the low complexity akin to the PWM, while also retaining the high accuracy of the SWM, thereby conferring significant advantages.

3.1 Closed-form expression of SNR↩︎

For downlink communication, considering the free-space LoS propagation, the user’s received signal can be expressed as

where \(P\) and \(s\) are the transmit power and information-bearing signal of the user, respectively; \(\boldsymbol{\eta} \sim \mathcal{CN}(0, \sigma^2 \mathbf{I}_{MN})\) is the AWGN with covariance matrix \(\sigma^2 \mathbf{I}_{MN}\); \(\boldsymbol{f} \in \mathbb{C}^{MN \times 1}\) is the transmit beamforming vector, with \(\|\boldsymbol{f}\| = 1\). Note that, although the uplink system model in Sections II, IV, and V employs a single-RF-chain analog combiner at the UAV receiver, we adopt a downlink model with digital beamforming here to facilitate tractable analysis. This allows us to derive closed-form SNR expressions and evaluate the upper-bound performance under ideal beamforming strategies. Thus the resulting received SNR can be given by \(\gamma = P \left| \boldsymbol{g}(\theta, \varphi, r) \boldsymbol{f} \right|^2/\sigma^2.\)

According to the maximal-ratio combining (MRC) strategy[[18]][34], i.e., \(\boldsymbol{f}^* = \frac{\boldsymbol{g}(\theta, \varphi, r)}{\|\boldsymbol{g}(\theta, \varphi, r)\|},\) the received SNR will reach its maximum value, which is

where \(\bar{P} = \frac{P}{\sigma^2}\) is the transmit SNR.

None

Following from (2) and (6 ), the SNR expression for SWM can be given by (4).

Proposition 1. The expression of the SNR for SWM can be obtained by

where

and \(S=A_{px}+\frac{D_x d}{2}+\frac{d}{2}\), which \(A_{px}\) represents the array aperture alone the \(x\)-axis, and \(\varrho\) is the argument of the function \(h(\cdot)\).

Proof: Please refer to Appendix A of the supplemental material.

Following from (3) and (6 ), the SNR expression for PWM can be calculated as

Incorporating (5 ), (?? ) and (?? ) into (6 ), the SNR expression for HSPWM can be given by (5).

None

To obtain a clearer analytical solution, the scenario is simplified by setting \(M_y=1\) and \(N_y=1\). Consequently, \(M\) UAVs are arranged in a row alone the \(x\)-axis, each carrying a ULA with \(N=N_x \times 1\) elements. Under these conditions, we derive the analytical SNR expression for three models.

Proposition 2. The analytical SNR expression for SWM under the case where the UAV carries a ULA can be derived as

Proof: Please refer to Appendix B of the supplemental material.

Remark: It is worth noting that (9 ) is structurally similar to (7 ). The difference lies in the fact that the ULA is less significantly affected by the distance \(r\) in the context of SWM compared to \(\gamma_{\mathbf{SWM}}\) in (7 ). From (9 ), it can be observed that the maximum SNR for multi-UAV distributed near-field communication is determined by the geometric configuration of the array, such as the total array size \(S\) and the UAV spacing \(D_x d\), as well as the distance \(r\) between the user and \(U_1\), the azimuth angle \(\theta\) and elevation angle \(\varphi\).

The SNR expression for PWM when the UAV is equipped with a ULA remains consistent with (8 ), i.e., \(\gamma_{\mathbf{PWM,ULA}} = \frac{\bar{P} \beta_0 M_xN_x}{r^2 }\), which increases linearly with \(MN\) and becomes unbounded as \(MN\) approaches infinity. This result is clearly impractical, fundamentally attributable to the inapplicability of the PWM assumption when \(MN\) is large, leading the communication scenario into the near-field domain. Additionally, it is crucial to recognize that while \(\gamma_{\mathbf{SWM,ULA}}\) depends on both the distance \(r\) and the geometric configuration of the array, \(\gamma_{\mathbf{PWM,ULA}}\) is solely dependent on the distance \(r\).

Proposition 3. The analytical SNR expression for HSPWM can be calculated as

where \(S_D=(M_x-1)D_x d\).

Proof: Please refer to Appendix C of the supplemental material.

Remark: The \(S_D\) can be understood as the aperture of a \(M_x \times 1\) ULA with a spacing of \(D_xd\). Comparing (9 ) and (10 ), it is evident that, different from SNR which is based on the SWM, the SNR based on the HSPWM assumption is solely related to the aperture \(S_D\), and the impact of \(N\) on SNR is linear. In other words, by jointly comparing (8 ), (9 ) and (10 ), the HPSWM integrates the characteristics of both the PWM and the SWM. Not only does it exhibit a linear increase with the growth of \(N\), but also depends on the geometric configuration of the ELAA, such as the aperture \(S_D\), the UAV spacing \(D_xd\), as well as the distance \(r\) between the user and \(U_1\), the elevation and azimuth angles.

3.2 Performance analysis↩︎

To further investigate the superiority of HSPWM in the multi-UAV near-field system, numerical simulations are developed to evaluate the analytical results derived in the previous section.

To enhance the generality of the scenario, we have configured the UAVs with a 4\(\times\)​4 structure, i.e., \(M_x=M_y=4\). The number of antennas for each UAV is set to \(N_x=N_y=12\). The carrier frequency is 10 GHz, and the inter-antenna spacing within the UAV is \(d=\frac{\lambda}{2}=0.015\) m. The spacing between the UAVs alone the \(x\) and \(y\) directions is \(D_x d=D_y d=50d=0.75\) m, and the coordinates of the user is set as \(\boldsymbol{p}_E=[0, 0, 0]^T\).

Fig. 2 illustrates the SNRs \(\gamma_{\mathbf{PWM}}, \gamma_{\mathbf{SWM}}\) and \(\gamma_{\mathbf{HPSWM}}\) versus the UAV spacing \(D_x d\), with different angles \(\varphi\) and \(\theta\). As shown in Fig. 2, \(\gamma_{\mathbf{PWM}}\) is independent of \(D\), \(\varphi\) and \(\theta\), which is in accordance with the theoretical analysis. In contrast, the SNRs \(\gamma_{\mathbf{SWM}}\) and \(\gamma_{\mathbf{HPSWM}}\) exhibit different trends for various UAV orientations \(\varphi\) and \(\theta\). It is also observed that \(\gamma_{\mathbf{PWM}}\) surpasses \(\gamma_{\mathbf{SWM}}\) when \(\varphi = -30^{\circ}\) and \(\theta = 0^{\circ}\), but falls below \(\gamma_{\mathbf{SWM}}\) when \(\varphi = 45^{\circ}\) and \(\theta = 30^{\circ}\). This demonstrates the necessity for near-field analysis in ELAA communication. Additionally, it is noted that \(\gamma_{\mathbf{HSPWM}}\) closely aligns with \(\gamma_{\mathbf{SWM}}\) regardless of the spacing variation. According to the previous theoretical analysis, the HSPWM exhibits the lower theoretical complexity compared to the SWM and better reflects the SNR trend in the near field, offering significant advantages. Furthermore, as observed in Fig. 2 and indicated by (9 ), the \(\gamma_{\mathbf{SWM}}\) can be either higher or lower than its PWM counterpart, depending on the directional cosine \(\Psi\), which reflects the angular sensitivity of SWM.

Fig. 3 illustrates the SNRs \(\gamma_{\mathbf{PWM}}, \gamma_{\mathbf{SWM}}\) and \(\gamma_{\mathbf{HPSWM}}\) versus the element number \(M N\) by increasing the number of UAVs \(M\) with different angles \(\varphi\) and \(\theta\). As shown in Fig. 3, under the same conditions \(\varphi\) and \(\theta\), for the relatively small element number \(MN\), the three models increase with the same growth trend. However, as \(MN\) further increases, the SNR based on the PWM and the SNRs based on SWM and HSPWM exhibit two different growth trends. Furthermore, it is observed that the SNRs based on the SWM and HSPWM models show different trends for different UAV orientations \(\varphi\) and \(\theta\), while \(\gamma_{\mathbf{PWM}}\) is independent of \(\varphi\) and \(\theta\). Importantly, regardless of the element number variation, \(\gamma_{\mathbf{HPSWM}}\) coincides exactly with \(\gamma_{\mathbf{SWM}}\). Both Figs. 2 and 3 show that as the UAV spacing and the number of antenna elements change, the SNR gap between the PWM and the SWM/HSPWM increases, which confirms the necessity of near-field analysis in distributed ELAA. It is also worth noting that the HSPWM and SWM curves always coincide under practical UAV configurations³. This is because the number of antennas mounted on a single UAV is physically limited (e.g., a typical UAV frame size of 47 cm × 58 . 5 cm × 21 . 5 cm (e.g., based on DJI Matrice 30) only supports up to around 31 elements), making the planar-wave assumption within each subarray in HSPWM sufficiently accurate. Therefore, HSPWM not only maintains modeling precision but also offers reduced complexity, validating its suitability as an efficient substitute for SWM in practical distributed ELAA scenarios.

4 Channel Estimation↩︎

In this section, leveraging the inherent characteristics of the distributed near-field scenario of multi-UAV, we design effective near-field channel estimation methods utilizing the accurate SWM and its approximation HSPWM, i.e., SD-OMP and tensor-OMP. In our approach, we assume that each UAV receives the user’s pilot signal using its onboard analog combiner, and the resulting low-dimensional compressed signals are subsequently shared with a central fusion node (e.g., a leader UAV or ground station) for joint processing. This centralized assumption enables high-accuracy channel estimation and simplifies algorithm design. Investigation of distributed estimation schemes and the effects of practical inter-UAV communication constraints will be considered in future work.

4.1 SD-OMP based on SWM↩︎

In this section, we propose a sparse estimation method based on the SWM. This method extends the angle-distance estimation framework from the PD to the SD, enabling the estimation of elevation, azimuth, and distance, and thus generating the SD dictionary. Moreover, compressive sensing techniques are employed to extract the essential angle and distance parameters necessary for channel reconstruction.

Specifically, the dictionary representation of channel can be formulated as

where

\(\left\{\theta_1,\ldots,\theta_A\right\}\), \(\left\{\varphi_1,\ldots,\varphi_B\right\}\) and \(\left\{r_1,\ldots, r_C\right\}\) denote the sampling points for the azimuth angle \(\theta\), the elevation angle \(\varphi\) and the distance \(r\), respectively. \(\boldsymbol{G}\) denotes the SD dictionary, and \(U=ABC\) represents the number of atoms. \(\boldsymbol{z}\in \mathbb{C}^{U \times 1 }\) is a sparse vector, which physically signifies the reconstruction of \(\boldsymbol{h}\) using the minimal number of non-zero values. The non-zero values represent the channel path gains, and these non-zero value indices indicate the selected atoms for the channel reconstruction.

The problem is transformed into minimizing the number of non-zero entries of \(\boldsymbol{z}\), hence we have

where \(\epsilon\) denotes the precise factor.

The problem proposed in (11 ) can be addressed by integrating it with the traditional OMP algorithm for sparse signal recovery, namely SD-OMP-ongrid. Due to its universality in many literatures like [25], a detailed explanation has been omitted.

After recovering \(\boldsymbol{z}\), the restored channel \(\hat{\boldsymbol{h}}\) can be obtained according to (?? ). The channel parameters of all paths can be obtained by \((\widehat{\boldsymbol{\theta}},\widehat{\boldsymbol{\varphi}},\widehat{\boldsymbol{r}}) = [(\widehat{\theta}_1, \widehat{\varphi}_1, \widehat{r}_1), \ldots, (\widehat{\theta}_L, \widehat{\varphi}_L,\widehat{r}_L)]\).

Furthermore, to achieve more accurate channel parameters, a new optimization problem can be formulated as

where

is the estimation result by SD-OMP-ongrid.

The Nelder-Mead algorithm [35], a derivative-free method, is utilized to solve (12 ) by taking the SD-OMP-ongrid results as initial values. This approach, termed SD-OMP-offgrid, enables lightweight local refinement in low-dimensional angular space without requiring gradient information. It is particularly suitable for our non-linear near-field model with array perturbations, where analytical gradients are difficult to obtain.

4.2 Tensor-OMP based on HSPWM↩︎

The methods based on SD-OMP provide an effective solution for multi-UAV near-field channel estimation. However, due to their high computational complexity required by the large-scale dictionary, these methods do not fully leverage the characteristics of multi-UAV near-field scenario, thereby not suitable for distributed ELAA systems. As a result, we propose a sparse estimation method based on the HSPWM in this subsection. Leveraging the inherent properties of the HSPWM, the near-field array response can be represented as the outer product of three vectors. It expands the PD dictionary to encompass the tensor dictionary, capitalizing on its distinctive mathematical formulation to simplify the representation of the dictionary.

4.2.1 Tensor dictionary↩︎

Inspired by [36], the construction of the tensor atoms is necessary. Reviewing (?? ), it can be readily derived that

where \(\boldsymbol{u}(\theta, \varphi) = \left[1, \ldots, e^{j \frac{2\pi}{\lambda} d(N_x-1) sin\varphi sin\theta }\right]\in \mathbb{C}^{N_x \times 1 }\), \(\boldsymbol{v}(\theta, \varphi) = \left[1, \ldots, e^{j \frac{2\pi}{\lambda} d(N_y-1) sin\varphi cos\theta }\right]\in \mathbb{C}^{N_y \times 1 }\).

Combining (5 ) and (13 ), it is evident that \(\boldsymbol{g}_{\mathbf{HSPWM}}(\theta, \varphi, r)\) can be represented by the Kronecker product of three vectors, which can be reshaped into a tensor with a closed-form expression. According to the rule of tensor vectorization, it is given by

where "\(\circ\)" denote the outer product, \(\boldsymbol{\mathcal{G}}(\theta, \varphi, r)\) is a third-order tensor with the size of \(N_y \times N_x \times M\) and devec(·) denotes the devectorization operator.

Consequently, the tensor-form channel can be expressed as:

The dictionary representation of tensor-form channel reconstruction can be formulated as:

where \(\boldsymbol{\mathcal{T}}\in \mathbb{C}^{A\times B \times C \times N_y \times N_x \times M }\) is a sixth-order tensor, which constitutes the tensor dictionary. \(\boldsymbol{\mathcal{T}}_{a,b,c}=\boldsymbol{\mathcal{G}}(\theta_{a},\varphi_{b}, r_{c})=\boldsymbol{v}(\theta_a, \varphi_b) \circ \boldsymbol{u}(\theta_a, \varphi_b) \circ \boldsymbol{b}(\theta_a, \varphi_b, r_c) \in \mathbb{C}^{N_y \times N_x \times M }\) is the \((a,b,c)\)-th atom, where \(\theta_a\), \(\varphi_b\) and \(r_c\) denote the \(a\)-th azimuth angle sample, the \(b\)-th elevation angle sample and the \(c\)-th distance sample of the tensor dictionary. \(\boldsymbol{\kappa}\) is used to store \(\kappa_{a,b,c}\), which is a sparse vector. If the value of \(\kappa_{a,b,c}\) is non-zero, the \((a,b,c)\)-th atom from the tensor dictionary is selected and the corresponding value represents the channel gain.

According to (3 ), the received signal can be rewritten as

where \(\overline{\boldsymbol{u}}(\theta_l, \varphi_l)=\boldsymbol{W}_x^H \boldsymbol{u}(\theta_l, \varphi_l)\in \mathbb{C}^{Q_x \times 1}\), \(\overline{\boldsymbol{v}}(\theta_l, \varphi_l)=\boldsymbol{W}_y^H \boldsymbol{v}(\theta_l, \varphi_l)\in \mathbb{C}^{Q_y \times 1}\). The result is obtained using \((\mathbf{A} \otimes \mathbf{B})(\mathbf{C} \otimes \mathbf{D}) = \mathbf{AC} \otimes \mathbf{BD}\). Additionally, vec(·) denotes the vectorization operator.

Employing the devectorization operation, the tensor-form received signal can be expressed as

where \(\boldsymbol{\mathcal{Y}}\in \mathbb{C}^{Q_y\times Q_x \times M }\) is the third-order tensor.

According to (14 ) and (15 ), we define \(\overline{\boldsymbol{\mathcal{WT}}}_{a,b,c}=\overline{\boldsymbol{v}}(\theta_a, \varphi_b) \circ \overline{\boldsymbol{u}}(\theta_a, \varphi_b) \circ \boldsymbol{b}(\theta_a, \varphi_b, r_c)\). Therefore, the dictionary representation of the received signal can be expressed as

The channel reconstruction problem is then transformed into a tensor recovery issue utilizing a tensor dictionary, which is given by

After recovering \(\boldsymbol{\kappa}\), the restored channel \(\hat{\boldsymbol{\mathcal{H}}}\) can be obtained according to (14 ).

Since (?? ) requires computation with a tensor dictionary, its recovery process deviates from the standard method outlined in (11 ). In this subsection, we present a greedy solution that includes the selection of the tensor atom to tackle this distinctive challenge.

4.2.2 3D Projection↩︎

The projection of the received signal \(\boldsymbol{\mathcal{Y}}\) onto the third-order tensor is defined as:

where

with \(\times_1, \times_2, \times_3\) denoting the \(n\)-mode product of a tensor with the vector.

A deeper understanding of the tensor multiplication process can be achieved by Fig. 4, which a grid represents a third-order tensor. The candidate atoms of Fig. 4 are extracted and reassembled into \(\overline{\boldsymbol{\mathcal{WT}}}\), with dimensions of length \(A\), width \(B\), and height \(C\), thus forming a sixth-order tensor. A particular element within \(\overline{\boldsymbol{\mathcal{WT}}}\), denoted as \(\overline{\boldsymbol{\mathcal{WT}}}_{a,b,c}\) , has dimensions of \(Q_y \times Q_x \times M\). The multiplication of \(\boldsymbol{\mathcal{Y}}\) with \(\overline{\boldsymbol{\mathcal{WT}}}_{a,b,c}\) can be regarded as the sequential mode products of the tensor with three vectors. Each mode product reduces the order of the tensor by one, ultimately transforming it into a zeroth-order tensor, which is a scalar. The final result consists of \(ABC\) values. Identifying the indices of the maximum value among these values corresponds to the angles and distances that we are seeking to determine.

4.2.3 Least squares↩︎

The set of selected atom indices is defined as \(\Gamma\). These selected atoms can be used to construct an approximation of \(\boldsymbol{\mathcal{Y}}\), that is:

We define a tensor \(\boldsymbol{\mathcal{R}}\in \mathbb{C}^{Q_y \times Q_x \times M}\) to represent the residual: \(\boldsymbol{\mathcal{R}}=\boldsymbol{\mathcal{Y}}-\boldsymbol{\mathcal{Y}}_{\Gamma}\). By seeking the optimal \(\boldsymbol{\kappa}_{\Gamma}\) to minimize the Frobenius norm of \(\boldsymbol{\mathcal{R}}\), the problem is essentially equivalent to a least squares problem. Consequently, the optimization problem can be formulated as

Unlike the traditional residual computation method used in OMP , we need to solve for the tensor-form coefficients, which undoubtedly poses a challenge to the solution process. This problem can be simplified through the operation of vectorization. Firstly, we have

where \(\boldsymbol{X} \in \mathbb{C}^{MQ \times length(\Gamma)}\) is constructed by the Kronecker product of three vectors, and \(length(\Gamma)\) denotes the number of selected atom.

Then the tensor-form minimization problem can be simplified to

\(\boldsymbol{\kappa}_{\Gamma}\) can be calculated as

After \(\boldsymbol{\kappa}_{\Gamma}\) is obtained, the residual \(\boldsymbol{\mathcal{R}}\) can be obtained as

4.2.4 Offgrid↩︎

The aforementioned procedures have been summarized into Algorithm 1, which is called tensor-OMP-ongrid.

Building upon these foundations, the pursuit of enhanced channel estimation precision leads us to further optimize the method. This advancement is realized through the tensor-OMP-offgrid, which evolves from the outcomes of the tensor-OMP-ongrid algorithm. Specifically, a new optimization problem can be formulated as

where

is the estimation result by the tensor-OMP-ongrid.

The Nelder-Mead algorithm is utilized to solve equation (?? ) by taking the results of tensor-OMP-ongrid as initial values. The detailed procedural steps for this refined approach are delineated comprehensively within Algorithm 1.

4.3 Disscussion↩︎

4.3.1 Sparsity discussion↩︎

(?? ) is undoubtedly a sparse recovery problem, which can be understood through Fig. 5.

In Fig. 5, each grid represents a third-order tensor, which corresponds to a atom. We first construct \(\overline{\boldsymbol{U}} = \left[ \overline{\boldsymbol{u}}(\theta_1, \varphi_1), \ldots, \overline{\boldsymbol{u}}(\theta_A, \varphi_B) \right] \in \mathbb{C}^{N_x \times AB}, \overline{\boldsymbol{V}} = \left[ \overline{\boldsymbol{v}}(\theta_1, \varphi_1), \ldots, \overline{\boldsymbol{v}}(\theta_A, \varphi_B) \right] \in \mathbb{C}^{N_y \times AB}, \boldsymbol{B} = \left[ \boldsymbol{b}(\theta_1, \varphi_1, r_1), \ldots, \boldsymbol{b}(\theta_A, \varphi_B, r_C) \right] \in \mathbb{C}^{M \times ABC}\), from which it can be inferred that the size of the cubic structure is \(AB\times AB\times ABC\). This also implies that the sampling required would be \(A^3B^3C\), which is substantially excessive. However, it is evident that due to the consistency of \(\theta\) and \(\varphi\) across the three vectors \(\overline{\boldsymbol{u}}\), \(\overline{\boldsymbol{v}}\) and \(\boldsymbol{b}\), the candidate atoms will only appear in specific positions, denoted by \(s=x=a\), \(t=y=b\). These positions are represented by the yellow grids in Fig. 5, which in this context signify \(\overline{\boldsymbol{\mathcal{WT}}}_{a,b,c}\). Notably, these yellow grids lie along the body diagonal of the tensor cube, ensuring the angular consistency \((\theta_a, \varphi_b)\) across the three components \(\overline{\boldsymbol{u}}\), \(\overline{\boldsymbol{v}}\), and \(\boldsymbol{b}\). The ellipsis “\(\cdots\)” along the third dimension denotes the continuation over the range index \(r_c\), showing that only atoms with angular alignment contribute to the sparse representation. The remaining blue grids indicate invalid or impossible atoms that violate the angular consistency rule, e.g., \(\overline{\boldsymbol{v}}(\theta_1, \varphi_4) \circ \overline{\boldsymbol{u}}(\theta_8, \varphi_2) \circ \boldsymbol{b}(\theta_3, \varphi_7, r_1),\) which combine unmatched angles and hence are excluded from the valid dictionary. This design effectively avoids redundant angular sampling and highlights the inherent sparsity of the problem.

4.3.2 Complexcity analysis↩︎

In SD-OMP, the projection step requires multiplying the large matrix \(\boldsymbol{W}^H \boldsymbol{G} \in \mathbb{C}^{MQ \times U}\) and computing the inner product with the residual vector. The projection step requires computing the matrix product \(\boldsymbol{W}^H \boldsymbol{G} \in \mathbb{C}^{MQ \times ABC}\), resulting in a complexity of \(\mathcal{O}(NM\cdot MQ \cdot U)\). In contrast, Tensor-OMP avoids explicit dictionary construction and performs three mode-\(n\) tensor products per atom, with each projection costing only \(\mathcal{O}(Q_x Q_y M)\). Therefore, the total projection complexity is reduced to \(\mathcal{O}( M Q\cdot U)\). This approach leverages the Kronecker structure and compressive projections, significantly reducing computational cost. While existing OMP variants such as SDAR-OMP [37] and 2D-OMP [38] have been proposed to reduce the complexity of sparse recovery, they differ significantly from our tensor-OMP method. Specifically, successive decision-aided recovery (SDAR)-OMP recovers multi-dimensional sparse vectors in a sequential manner under the Kronecker sensing model, which requires intermediate decision steps and lacks global tensor structure exploitation. 2D-OMP, on the other hand, operates in the matrix domain and is limited to 2D signal recovery. In contrast, our method leverages the structure of the near-field HSPWM array to define tensor atoms, and performs matching directly in the tensor domain, enabling batch projection onto the entire dictionary with significantly lower complexity.

5 Simulation results↩︎

In this section, several numerical simulations are conducted to demonstrate the effectiveness of the proposed method under various system parameters. The system operates at a central frequency of 10 GHz. The multi-UAV configuration is \(M_x=4\), \(M_y=2\), with the total number of RF chains being \(M=M_x M_y=8\). Each UAV is equipped with UPA of \(N_x=N_y=12\). The inter-antenna spacing within the UAV is \(d=\lambda/2\), and the spacing between the UAVs in the \(x\) and \(y\) directions is \(D_xd=D_yd=50d\). Consequently, the system aperture \(A_p\) is calculated as \(A_p=\sqrt{A_{px}^2 + A_{py}^2}\), where \(A_{py}=(M_y-1) D_y d+ (N_y-1) d\) denotes the aperture alone the \(y\) axis. So the Rayleigh distance can be computed as \(R=2 A_p^2/\lambda=444.63\) m. Furthermore, the azimuth angle \(\cos \theta\), the elevation angle \(\sin \varphi\) and the distance \(r\) follow the uniform distribution \(\mathcal{U}\left(-0.75, 0.75\right)\), \(\mathcal{U}\left(-0.75, 0.75\right)\) and \(\mathcal{U}\left(50, 80\right)\), respectively. Meanwhile, when constructing the dictionary, the range of the sampled azimuth angle \(\sin\varphi\) is \([-0.75, 0.75]\), the range of the sampled elevation angle \(\cos\theta\) is \([-0.75, 0.75]\), with the grid resolution set at \(2/N\), and the range of sampled distance \(r\) is [50 m, 80 m]. The distance sampling method in the polar-domain dictionary construction follows the non-uniform strategy proposed in [25], where the sampling is uniformly performed on \(1/r\) instead of \(r\), to reduce the column coherence and improve near-field sparsity, we have

The step size \(\Delta\) is set to 0.0003. ⁴

The benchmark scheme used in the study is described as follows:

• SD-OMP-far-field: The algorithm can be understood as the application of the SD-OMP-offgrid in the far field context. \(r\) is set as \(r=10^5 \text{m} \gg R\).

• SD-OMP-ongrid: As shown in .A.

• SD-OMP-offgrid: As shown in IV.A. Integrating the offgrid components in Algorithm 1 based on the SD-OMP-ongrid algorithm.

• Tensor-OMP-ongrid: As shown in Algorithm 1. \(\widehat{\boldsymbol{h}}\) is the channel reconstructed by tensor-OMP-ongrid.

• Tensor-OMP-offgrid: As shown in Algorithm 1. \(\widetilde{\widehat{\boldsymbol{h}}}\) is the channel reconstructed by tensor-OMP-offgrid.

• LS: Using the least squares (LS) estimator as an upper bound.

As depicted in Fig. 6, we present the normalized mean square error (NMSE) for various methods under conditions where \(L = 1\), \(Q_x = Q_y \in \{6, 8, 10, 12\}\), and the \(\text{SNR} \in \{-10, -5, 0, 5, 10, 15, 20\}\) dB. NMSE can be calculated as \(\text{NMSE} = \frac{\mathbb{E}(\| \boldsymbol{h} - \hat{\boldsymbol{h}}\|_2^2)}{\mathbb{E}(\| \boldsymbol{h} \|_2^2)}\). It is evident from Fig. 6(a) that at lower SNR levels, both tensor-OMP-ongrid and tensor-OMP-offgrid exhibit slightly higher NMSE compared to SD-OMP-ongrid and SD-OMP-offgrid. However, as the SNR increases, SD-OMP-offgrid and tensor-OMP-offgrid demonstrate superior performance over their ongrid counterparts. Notably, SD-OMP-ongrid shows the poorest performance, while tensor-OMP-offgrid achieves the best result. It is important to observe that with the increase of SNR, tensor-OMP-ongrid approaches a performance similar to that of SD-OMP-ongrid, highlighting the superiority of the tensor-OMP algorithm proposed in this paper. In Fig. 6(b), similar conclusions can be drawn. In Fig. 6(c) and Fig. 6(d), with the increase in the number of measurements \(Q_x = 10, 12\), SD-OMP-offgrid shows better performance. Although tensor-OMP-ongrid exhibits inferior performance, comparatively, regardless of the variation in \(Q_x\), tensor-OMP-offgrid consistently outperforms tensor-OMP-ongrid in reducing NMSE more significantly than SD-OMP-offgrid does compared to SD-OMP-ongrid. This suggests that tensor-OMP-offgrid is more effective in minimizing NMSE across different values of \(Q\).

In Fig. 7, we plot the NMSE for various methods when \(L=2\). At this point, considering the characteristics of UAV communication, channel model (4 ) comprises one LOS path and one Non-LOS (NLoS) path. As the number of paths \(L\) increases, the complexity of the recovery problem also increases. The enhanced sparsity leads to a greater number of parameters that need to be estimated, which in turn causes changes in the NMSE for all methods. In general, Fig. 7 shows similar performance to that of Fig. 6. It can be observed that SD-OMP-ongrid performs better than tensor-OMP-ongrid when \(L=2\). On this basis, SD-OMP-offgrid achieves a greater improvement in performance, reflecting the superiority of the SD-OMP algorithm. However, as the SNR increases, the tensor-OMP-offgrid algorithm, in most cases, realizes performance comparable to that of SD-OMP-offgrid, making it suitable for practical applications. Notably, as SNR increases, the NMSE curve gradually flattens, due to the limited number of measurements, which prevents NMSE from approaching zero even as SNR tends to infinity.

6 Conclusion↩︎

In this paper, we investigated the SNR analysis and channel estimation for multi-UAV near-field systems, which can be regarded as a distributed ELAA architecture. By introducing the HSPWM model, we demonstrated through SNR analysis that it offers a simplified yet accurate framework for near-field channel modeling. For channel estimation, we developed SD-OMP and tensor-OMP approaches. SD-OMP extends traditional PD-OMP from angle-distance to elevation-azimuth-distance estimation, resulting in both ongrid and offgrid variants. Based on HSPWM, we further constructed a tensor dictionary and proposed tensor-OMP-on/offgrid algorithms. Simulation results show that tensor-OMP-offgrid achieves a good trade-off between computational complexity and estimation accuracy. As UAV jitter and synchronization errors are inevitable in practical deployments, understanding and mitigating their impact on system performance will be the focus of our future research.

7 Proof of (7 )↩︎

We first define a continuous function \(f(x,y,s,t)\) about \(x,y,s,t\) over the 4-directional hypercube, where \(x \in \left[-\frac{\xi}{2}, (M_x-\frac{1}{2}) \xi \right]\), \(y \in \left[-\frac{\xi}{2}, (M_y-\frac{1}{2}) \xi \right]\), \(s \in \left[-\frac{\xi}{2}, (N_x-\frac{1}{2}) \xi \right]\) and \(t \in \left[-\frac{\xi}{2}, (N_y-\frac{1}{2}) \xi \right]\), shown at (13). Given that \(\xi \ll 1\), within each unit hypercube, where \(\forall x \in \left[ \left( m_x - \frac{1}{2} \right) \xi, \left( m_x + \frac{1}{2} \right) \xi \right], y \in \left[ \left( m_y - \frac{1}{2} \right) \xi, \left( m_y + \frac{1}{2} \right) \xi \right], s \in \left[ \left( n_x - \frac{1}{2} \right) \xi, \left( n_x + \frac{1}{2} \right) \xi \right], t \in \left[ \left( n_y - \frac{1}{2} \right) \xi, \left( n_y + \frac{1}{2} \right) \xi \right]\), we have \(f(x,y,s,t)\approx f(m_x\xi,m_y\xi,n_x\xi,n_y\xi)\).

None

Based on the definition of Riemann sums, we have

By substituting (13) into (16 ), we have (14), shown at the top of the page, where (a) holds according to \(\int \frac{dx}{A + 2Bx + Cx^2} = \frac{1}{\sqrt{AC - B^2}} \arctan\left(\frac{Cx + B}{\sqrt{Ac - B^2}}\right)\) and (b) holds from the integral formulas \(\int \arctan\left(\frac{x}{a}\right) dx = x \arctan\left(\frac{x}{a}\right) - \frac{a}{2} \ln\left(a^2 + x^2\right)\) in [39], respectively. From (14) and (4), we have (7 ). The proof of (7 ) is completed.

None

8 Proof of (9 )↩︎

We first give the analytical expression of \(\gamma_{\mathbf{SWM,ULA}}\) under the condition that the UAV carries ULA.

Similar with Appendix A, according to the principles of (a) and (b), we have (15).

None

9 Proof of the (10 )↩︎

The SNR expression of \(\gamma_{\mathbf{HSPW,ULA}}\) under the condition that UAV carries ULA can be given by

The remaining proof of (10 ) is similar to (a) of Appendix A, which is omitted for brevity.

References↩︎

---

1. Tianyu Huo is with the Institute of Image Communication and Network Engineering, Shanghai Jiao Tong University, Shanghai 200240, China. (e-mail: huotianyu@sjtu.edu.cn).

   Jian Xiong is with the Institute of Image Communication and Network Engineering, Shanghai Jiao Tong University, Shanghai 200240, China; he is also with Zhengzhou Industrial Technology Research Institute of SJTU (e-mail: xjarrow@sjtu.edu.cn).

   Yiyan Wu is with the Department of Electrical and Computer Engineering, Western University, London, Ontario, Canada. (email: Yiyan.wu@ieee.org)

   Songjie Yang is with the National Key Laboratory of Wireless Communications, University of Electronic Science and Technology of China, Chengdu 611731, China. (e-mail: yangsongjie@std.uestc.edu.cn).

   Bo Liu is with the School of Computer Science, University of Technology Sydney, Ultimo, NSW 2007, Australia. (e-mail: bo.liu@uts.edu.au).

   Wenjun Zhang is with the Department of Electronic Engi- neering, Shanghai Jiao Tong University, Shanghai 200240, China. (e-mail:zhangwenjun@sjtu.edu.cn).↩︎

2. All UAVs are assumed to upload their received pilot signals to a central node (e.g., a ground station or a leading UAV) for centralized processing[30]. The centralized scheme assumes that UAVs only upload a small set of sparse channel parameters (e.g., user position and path gains) during predefined time slots, and the central node reconstructs the channel accordingly [31].↩︎

3. The overlap of HSPWM and SWM curves results from the limited number of antennas per UAV, which ensures the far-field assumption within each subarray remains valid in practice.↩︎

4. Distance Sampling: The chosen step size satisfies the condition \(\Delta < \frac{1}{Z_\Delta}\) to control dictionary coherence, where \(Z_\Delta = D^2/2 \beta^2 \lambda.\) For our system setup with estimated \(Z_\Delta \approx 1235\)( the desired coherence leve \(\beta=0.03\)), the maximum allowable step size is \(\Delta_{\text{max}} \approx 0.00081\), validating our use of \(\Delta = 0.0003\).↩︎