Energy and Information Beamforming in Airborne Massive MIMO System for Wireless Powered Communications

Energy supply and information backhaul are critical problems for wireless sensor networks deployed in remote places with poor infrastructure. To deal with these problems, this paper proposes an airborne massive multiple-input multiple-output (MIMO) system for wireless energy transfer (WET) and information transmission. An air platform (AP) equipped with a two-dimensional rectangular antenna array is employed to broadcast energy and provide wireless access for ground sensors. By exploiting the statistical property of air-terrestrial MIMO channels, the energy and information beamformers are jointly designed to maximize the average received signal-to-interference-plus-noise ratio (SINR), which gives rise to a statistical max-SINR beamforming scheme. The scheme does not rely on the instantaneous channel state information, but still requires large numbers of RF chains at AP. To deal with this problem, a heuristic strongest-path energy and information beamforming scheme is proposed, which can be implemented in the analog-domain with low computational and hardware complexity. The analysis of the relation between the two schemes reveals that, with proper sensor scheduling, the strongest-path beamforming is equivalent to the statistical max-SINR beamforming when the number of AP antennas tends to infinity. Using the asymptotic approximation of average received SINR at AP, the system parameters, including transmit power, number of active antennas of AP and duration of WET phase, are optimized jointly to maximize the system energy efficiency. The simulation results demonstrate that the proposed schemes achieve a good tradeoff between system performance and complexity.


Motivations
In recent years, wireless energy transfer (WET) has been regarded as a promising technology to increase the battery-lifetime of energy-constrained wireless sensor nodes [1,2]. The WET technology allows sensor nodes to harvest energy from surrounding electromagnetic radiation instead of wired energy sources. Many works investigated the WET networks where the traditional ground base station (BS) or access point also acts as an energy station, which can charge the sensor nodes through WET [3]. However, such a configuration is still costly or even infeasible if the sensors are deployed in remote places with poor infrastructure (that may be used for forest fire detection, climate monitoring, etc.).

Focus and Contributions
In this paper, to deal with the energy supply and information backhaul problems for wireless sensor networks in remote places, we propose an airborne massive MIMO system for WET by combining the AP with massive MIMO technology. The system consists of a low-altitude AP equipped with a two-dimensional (2D) rectangular antenna array and a number of sensors on the ground, as shown in Figure 1. The transmission phase is divided into a WET phase during which the AP charges sensors through downlink WET and an uplink wireless information transmission (WIT) phase during which the AP collects data sent by sensors. In the considered hybrid air-terrestrial network, the biggest challenge is the design of practical energy and information beamforming schemes. Different from the BS of a terrestrial network, AP is usually hardware-and computational complexity-limited due to the consideration of weight, fabricating cost and energy supply [8]. Traditional linear beamforming schemes, such as match filter (MF) beamforming [10] and zero-forcing beamforming [12], require a large number of radio frequency (RF) chains at AP and have high computational complexity. This may make them infeasible when applied on AP. Another challenge lies in the acquisition of instantaneous channel state information (CSI) at the AP. Due to the mobility of AP, the sensor must transmit the uplink pilot sequence frequently, which causes high training overhead and energy consumption. This paper studies the three-dimensional (3D) energy and information beamforming scheme for an airborne massive MIMO system considering the above challenges. The main contributions are as follows. • By exploiting the statistical property of air-terrestrial MIMO channels, the energy and information beamformers are jointly designed to maximize the average received signal-to-interference-plus-noise ratio (SINR), which gives rise to a statistical max-SINR scheme.
Although the scheme does not rely on instantaneous CSI, its implementation still requires large numbers of RF chains at AP. When the number of RF chains is limited, a heuristic strongest-path energy and information beamforming scheme is proposed. The scheme requires only virtual angle of departure (AoD) information of sensors and can be implemented in the analog-domain with low hardware complexity. The analysis of the relation between two schemes reveals that, with proper sensor scheduling, the strongest-path beamforming is equivalent to the statistical max-SINR beamforming when the number of AP's antennas tends to infinity.

•
Based on the asymptotic approximation of average received SINR at AP, the system parameters, including transmit power, number of active antennas of AP and duration of the WET phase, are jointly optimized to maximize the system EE under the proposed strongest-path beamforming scheme. • Numerical simulations are presented to evaluate the spectral and energy efficiency performances of the proposed beamforming schemes under different system parameters. The results demonstrate the superiority of the proposed schemes in the airborne massive MIMO systems.

Related Literature
The modeling and application of hybrid air-terrestrial networks for access of ground users have drawn much research attention recently. A statistical propagation model has been introduced in [14] to predict the air-to-ground path loss between AP and ground users based on ray tracing simulation. The authors in [15] considered the modeling of small-scale fading for air-to-ground channels, where both line-of-sight (LoS) and non-LoS (NLoS) channel environments were considered. In [4], the authors discussed the basic networking architecture and major design challenges in UAV-aided hybrid air-terrestrial networks. In [16], the authors investigated the 3D beamforming design for a hybrid air-terrestrial network with a large-scale antenna array and analyzed the effect of AoD imperfection on the system performance.
On the other hand, the WET in terrestrial networks has also been extensively investigated for both indoor and outdoor scenarios [17][18][19][20][21][22][23][24][25][26][27]. In [17,18], the authors investigated the outage performance of the indoor WET system and showed that increasing the variance of the log-normal channel would degrade the system performance. In outdoor scenarios, to combat the path loss due to the long charging distance, a number of works considered the utilization of multi-antenna technologies in the WET system. In particular, the benefits and challenges of combining WET with multi-antenna technologies were discussed in [19,20]. The authors in [21] considered the EE optimization in a massive MIMO system with WET and proposed an iterative algorithm to compute the optimal transmit power and duration of the WET phase. The authors in [22] investigated the throughput optimization problem for a WET-enabled massive MIMO system. In [23], an asymptotically optimal downlink power allocation strategy was proposed to maximize the uplink sum rate in a WET-enabled massive MIMO system with linear beamforming. In [24], the 3D massive MIMO was utilized for WET, where the BS applied linear beamforming for energy broadcasting and the users applied power splitting for information detection and energy harvesting. In [25], the energy and information beamforming design for a full-duplex massive MIMO system was investigated. The authors in [26,27] investigated the combined benefits of WET and cooperative relaying systems with large-scale antenna arrays. Due to practical constraints, [19][20][21][22][23][24][25][26][27] have assumed that the energy harvesters (e.g., sensors) are equipped with a single antenna. Recently, several works began to drop this assumption and consider multiple-antenna sensors [28][29][30]. It was shown that, by equipping more antennas at the sensor, the energy harvesting and achievable rate performances can be improved through coherent processing at the sensor side. However, the increase of antennas also incurs problems such as high circuit energy consumption and hardware complexity, which may be unaffordable for sensors. Therefore, in this paper, we assume a single-antenna sensor, and the general case with multiple-antenna sensors will be considered as future work.

Organization and Notations
The rest of the paper is organized as follows. Section 2 presents the system model. Section 3 presents the energy and information beamforming design. Section 4 presents the EE optimization scheme. Section 5 presents the simulation results. Section 6 draws the conclusions.

Channel Model
Let h l ∈ C N x N y ×1 denote the channel vector between AP and sensor l. Using the air-terrestrial 3D MIMO channel model in [15,31], h l can be expressed as: (1) In (1), h LoS l ∈ C N x N y ×1 denotes the deterministic LoS channel component and h NLoS l ∈ C N x N y ×1 is the random vector characterizing the non-LoS channel component, which are given by: with: where β l and K l denote the large-scale fading and Rician factor, respectively. ϕ l ∈ [−π/2, π/2] and θ l ∈ (0, π/2] denote the AoDs of sensor l in the horizontal and vertical directions, and ∆ϕ l and ∆θ l denote the corresponding angular spreads. Since AP is much higher than ground obstacles, there is little scatter around it. In this case, the physical signal propagation paths are mainly caused by the scattering process in the vicinity of sensors. This will result in very narrow angular spreads ∆ϕ l and ∆θ l [11]. d x and d y denote the antenna spacings along the x axis and y axis, respectively. λ denotes the carrier wavelength, and r l (θ, ϕ) denotes the complex response gain associated with the direction (θ, ϕ). We assume that r l (θ, ϕ) has zero-mean, i.e., E[r l (θ, ϕ)] = 0, and the response gains for different angles are uncorrelated, i.e., where S l (ϕ, θ) represents the channel power angle spectrum (PAS), which characterizes the channel power distribution in the angular domain. Define ρ x = d x λ cos θ sin ϕ and ρ y = d y λ cos θ cos ϕ. Using a similar procedure as that in [16], the channel correlation matrix can be expressed as: where ρ x,l = d x λ cos θ l sin ϕ l and ρ y,l = d y λ cos θ l cos ϕ l denote the virtual AoDs corresponding to the sensors' physical AoDs in the horizontal and vertical directions. S l ρ x , ρ y denotes the PAS with respect to the virtual AoDs ρ x and ρ y . The integral boundaries are given by:

Energy Transfer and Information Transmission
We adopt frame-based transmission where the length of one frame is denoted by T. Each frame consists of a WET phase and a WIT phase.
In the WET phase, the AP charges sensors via downlink energy beamforming. The duration of WET phase is τT with 0 < τ < 1. According to the law of energy conservation, the harvested energy at sensor l can be expressed as [21][22][23] (As in [21][22][23], we have neglected the non-linearity in the energy harvesting process. This is a good assumption for long-range WET, as shown in [32].): where w E,j denotes the energy beamformer for sensor j. η (0 < η < 1) denotes the energy conversion efficiency, i.e., the ratio between the harvested energy and received energy. p j denotes the power allocated to the beamformer of sensor j, which satisfies ∑ L j=1 p j = p sum . The first term of the right-hand side indicates the expected energy that can be harvested at sensor l, and the second term indicates the energy leakage from the signals sent to other sensors.
Note that the schemes in [19,[21][22][23][24][25][26][27][28][29][30] assume that all harvested energy is used for the power amplifier of the sensor during the WIT phase. Different from the previous works, in this paper, we employ a more practical power consumption model for the sensor. In particular, after the WET phase, a part of the harvested energy is used for necessary processing at the sensor to achieve its main task, and the remaining energy is used for uplink transmission during the WIT phase. A similar model has also been applied in [20]. Let E p l denote the processing energy, in the WIT phase of time period (1 − τ)T, the transmit power of sensor l can be expressed as: Equation (9) has assumed that the harvested energy E l is greater than the processing energy E p l . Otherwise, P l should be set to zero, and an outage occurs.
To simplify the problem, we assume that the energy leakage, i.e., the second term of (8), contributes little to the total harvested energy. As will be seen later, this assumption is well satisfied under the proposed sensor scheduling scheme. Therefore, the transmit power of sensor l in (9) can be approximated as: The received signal at AP during the WIT phase is: where H = [h 1 , h 2 , · · · , h L ] and Λ = diag {P 1 , P 2 , · · · , P L }. x = [x 1 , x 2 , · · · , x L ] T , and x l denotes the unit-power transmit signal of sensor l. n denotes the additive white Gaussian noise (AWGN) vector with distribution CN 0, σ 2 I N x N y . To decode the signal of sensor l, AP combines the received signal by multiplying the information beamformer w D,l , i.e., For convenience, we normalize w E,l and w D,l as w E,l = w D,l = 1. The average received SINR of the transmission from sensor l to AP can be expressed (since the proposed scheme does not rely on the instantaneous CSI, we consider average received SINR as the performance metric): The approximation in the second line is based on the Mullens inequality [33], which is widely employed in the literature of massive MIMO. The approximation has been proven tight when the number of antennas (i.e., N x and N y ) is large [10].

Energy and Information Beamforming
In this section, we consider the energy and information beamforming design with statistical CSI. When large numbers of RF chains are available at AP, the energy and information beamformers are designed jointly to maximize the average received SINR at AP, which we call the statistical max-SINR beamforming. Although good performance is expected, the scheme causes high computational and hardware complexity when the number of AP's antennas becomes very large. To deal with this problem, we then propose a heuristic strongest-path energy and information beamforming scheme, which requires only virtual AoD information of sensors and can be implemented in the analog-domain with low hardware complexity. Moreover, with strongest-path beamforming, a simple analytical expression for average received SINR at AP, can be obtained. With this property, an efficient parameter optimization scheme can be developed to further improve the system performance, as will be seen in the next section. The analysis of the relation between two beamforming schemes reveals that, with proper sensor scheduling, the strongest-path beamforming is equivalent to the statistical max-SINR beamforming when the number of AP antennas tends to infinity.
In this section, for information beamforming design, we will assume that the harvested energy is enough, and as a result, the transmit power of sensor given by (10) is positive. When this assumption is not satisfied, the information beamforming problem becomes trivial since the received SINR at AP is always zero.

Statistical Max-SINR Scheme
With only statistical CSI, the beamforming problem is formulated as follows: Directly solving the above problem is difficult due to the lack of analytical expression of average received SINR. Moreover, the energy beamforming designs for different sensors are coupled through the multi-user interference (MUI) term (i.e., the first term on the denominator of (13)), which makes the problem more challenging.
In this subsection, instead of directly solving the average SINR maximization problem, we first propose a statistical signal-to-leakage-plus-noise ratio (SLNR) criterion to obtain a closed-form solution for the energy beamformer. As its name suggests, the idea is inspired by the conventional SLNR criterion based on instantaneous CSI (referred to as instantaneous SLNR), which has been proven near optimal for the downlink MIMO beamforming problem [34,35]. Then, by exploiting the asymptotic property of the channel covariance matrix, the information beamformer is obtained by solving the problem (14) optimally. Note that in the above strategy, w E,l and w D,l are solved jointly since we have not added any assumption to decouple the problem.

Energy Beamforming Design
Similar to the instantaneous SLNR [34,35], the statistical SLNR is defined as follows: where the first term on the denominator indicates the total power of interferences in the WIT phase caused by the transmission of sensor l. For convenience, we define G l,k = w H D,k h l h H l w D,k as the combining gain of information beamformer w D,k on the signal of sensor l. Substituting the expression of P l given by (10) into (15), the statistical SLNR can be rewritten as: where Based on (16), the energy beamformer is designed to maximize the lower bound of statistical SLNR, i.e., max It is easy to see that (17) is the generalized Rayleigh quotient problem, which can be solved as: Note that the matrices K l and have the same eigenvectors with C l . Let λ i denote the i-th largest eigenvalue of C l ; the i-th largest eigenvalue of K l can be expressed as: Therefore, (18) can be simplified as: As discussed in [22], the MF energy beamforming aims to maximize the harvested energy at the intended sensor with instantaneous CSI. It is interesting that the energy beamformer (20) under the statistical SLNR criterion gives rise to a similar interpretation in the sense of maximizing the average harvested energy at the intended sensor. Moreover, it is seen that the solution of the energy beamformer is decoupled with the information beamformer. This makes the closed-form solution for w D,l possible.

Information Beamforming Design
By substituting (20) into (14), the information beamforming problem can be expressed as follows: Again, (21) is in the form of generalized Rayleigh quotient problem, which can be solved as: Since the transmit power P l is correlated with the channel h l , there is no closed-form expressions for the expectations in (21). Therefore, the information beamforming based on (22) requires numerical evaluation with high computational complexity.
To overcome this challenge, below, we propose a closed-form solution for w D,l based on the asymptotic property of the sensor's transmit power P l . Lemma 1. When the number of AP antennas tends to infinity, i.e., N x , N y → ∞, the transmit power of sensor l converges to: Proof. See Appendix A.
Inserting (23) into (22), we get a closed-form solution for w D,l as follows: Even if closed-form expressions are available, the computation of beamformers in (20) and (24) has complexity O N 3 x N 3 y , which becomes challenging when the number of AP antennas grows large. Moreover, to implement the beamforming scheme, the required number of RF chains at AP is equal to N x N y . Since AP is hardware complexity limited, it is very important to reduce the number of required RF chains at AP for energy and information beamforming. This will be considered in the next subsection.

Strongest-Path Beamforming Scheme
In this subsection, we present a two-step heuristic approach to obtain a low-complexity beamforming scheme. In the first step, we design the energy and information beamformers jointly to maximize the average power of useful signal (i.e., the numerator of (13)), without considering the power of MUI (i.e., the first term on the denominator of (13)). In the second step, we propose an MUI-aware sensor scheduling scheme to mitigate the MUI under the beamforming scheme obtained in the first step.

Average Useful Signal Power Maximization
According to (13), the average power of useful signal can be expressed as: Directly maximizing Q U l with respect to w E,l and w D,l is still difficult. By using the inequality l can be obtained as: Since K l has the same eigenvectors as C l , designing the energy and information beamformers to maximize the lower bound gives rise to the solution as follows: By exploiting the asymptotic property of the channel covariance matrix, we can further approximate the energy and information beamformers using the virtual AoD information of the sensor. This is stated in the following theorem. Theorem 1. When the number of AP antennas tends to infinity, i.e., N x , N y → ∞, the energy and information beamformers converge to: and the lower bound of average useful signal power can be expressed as: Proof. This can be readily proven using the asymptotic expressions for the dominant eigenvector and eigenvalue of C l , presented in (A6) of Appendix A.
From Theorem 1, it is seen that the beamformers are designed to match the LoS path of the channel, i.e., to ensure that the largest beamforming gain is achieved at the LoS path. Moreover, for the commonly-used models, the PAS of NLoS channel S l ρ x , ρ y achieves the maximum at {ρ x , ρ y } = {ρ x,l , ρ y,l } [15,36], i.e., the virtual AoDs corresponding to the sensors' physical AoDs in horizontal and vertical directions. As a result, the proposed beamformers also match the strongest NLoS path in the angular domain. Therefore, this scheme is named the strongest-path energy and information beamforming. On the other hand, it is seen that both the energy and information beamformers have a constant envelope, which can be realized using phase shifting networks in the analog-domain. The total number of required RF chains at AP can be reduced to L, which is equal to the number of sensors [37].
Another important observation from Lemma 1 is that, if we neglect the processing energy E p l , the average useful signal power decays with the squared path loss. Thus, a low altitude of AP is preferred to improve the efficiencies of WET and WIT of the single sensor. However, with the decrease of altitude, the number of covered sensors by AP is reduced as well, which may degrade the performance of the whole hybrid air-terrestrial network. Therefore, the planning of AP's altitude is an interesting and important research topic. If the prior of sensor density is available, a raw idea to solve this problem is presented as follows. Note that for the given altitude of AP, the uplink achievable rate region of the sensor can be predicted based on (29) or the analytical SINR expression presented in Section 4. Moreover, the performance metric of the network, e.g., system spectral efficiency (SE), can be estimated using the sensor density information. With these in the hand, the altitude optimization can be formulated as a network performance maximization problem subject to the sensor quality of service (QoS) constraint, e.g., a threshold on the uplink achievable rate.

MUI Mitigation
According to (13) and Lemma 1, the power of MUI can be expressed asymptotically in the large {N x , N y } limit as: By exploiting the channel covariance matrix (5) and the expression of w D,l given by Theorem 1, Q MUI l can be rewritten as: A short derivation of (31) is presented in Appendix B. Before proposing the MUI-aware sensor scheduling scheme, we first present the following lemma, which determines the scaling behavior of MUI with respect to {N x , N y }. Based on Lemma 2 and the above discussion, we can build an MUI-aware sensor scheduling scheme as follows.

•
In the first step, we divide the whole region of virtual AoD along the x-direction into B x disjoint blocks, and the length of each block is much larger than max l ρ max x,l − ρ min x,l . Similarly, the whole region of virtual AoD along the y-direction is divided into the B y disjoint blocks, and the length of each block is much larger than max l ρ max y,l − ρ min y,l , as shown in Figure 2.
• In the second step, the sensors whose virtual AoDs ρ x,l , ρ y,l lie in the same small rectangular region in Figure 2 are gathered into one group. In this way, we will have at most B x B y sensor groups.

•
In the third step, for the low received SNR scenario, we pick one sensor from each group to serve a particular time-frequency resource. Thus, the maximum number of simultaneously served sensors is equal to B x B y . In the high received SINR scenario, the scheduled sensors are selected in the same way, but further divided into two clusters, as shown in Figure 2, where the sensors from red regions are assigned to Cluster 1 and those from white regions are assigned to Cluster 2. Each cluster is served on different time-frequency resources. Thus, the maximum number of simultaneously served sensors is equal to 1 2 B x B y . To ensure the fairness between all sensors, the remaining sensors in each group can be served on the other time-frequency resources using the same scheduling procedure.
Since in the first step, we have made the length of each block much larger than the spread of virtual AoD of each sensor, i.e., ρ max x,l − ρ min x,l and ρ max y,l − ρ min y,l , the proposed scheme can make Condition (a) and/or (b) satisfied approximately. In this case, the power of MUI is greatly mitigated. From the first and third steps, B x and B y can be viewed as the parameters that give a tradeoff between the allowed residual MUI and the number of simultaneously served sensors. Thus, considering the constraint in the first step, B x and B y should be designed to maximize some performance metric (e.g., the system SE). In practical application, since the feasible sets B x and B y are discrete and finite, this task can be completed by simple off-line searching algorithms based on the analytical results in Theorem 1 and (31).

Remark 1.
Note that both the strongest-path beamforming and sensor scheduling require the knowledge of virtual AoDs of sensors. By using the class of estimation methods based on discrete Fourier transform with zero padding [38], the virtual AoDs can be estimated at AP from the pilots transmitted by sensors with a certain quantification error.

Remark 2.
If the remaining battery power of a sensor is not enough to transmit the pilot signal, AP cannot obtain its AoD estimation. In this case, this sensor becomes a "dead sensor", which can no longer be severed.
To deal with this problem, the common signal of AP can be exploited. Usually, AP will broadcast periodically some common signals within its coverage region for the purpose of synchronization, access control, etc. The sensor with low battery power can harvest energy from these common signals until it has enough battery power to access the network actively. Some methods of designing an omnidirectionally common signal in a massive MIMO system can be found in [39].

Relation between Two Beamforming Schemes
According to (20) and (27) and Theorem 1, we can see that the energy beamformer of the statistical max-SINR scheme converges to that of the strongest-path beamforming in the large N x , N y limit. Moreover, if Conditions (a) and/or (b) are satisfied after the MUI-aware sensor scheduling, we can neglect the term incurred by MUI, i.e., ∑ L k=1,k =l (1−τ)T C k , in the information beamformer given by (24), and rewrite (24) as: Therefore, we can conclude that, under the MUI-aware sensor scheduling, the strongest-path beamforming is equivalent to the statistical max-SINR beamforming when the number of AP antennas tends to infinity.

Energy Efficiency Optimization
In this section, the transmit power, number of active antennas of AP and duration of the WET phase are jointly optimized to maximize the system EE. It is assumed that the low-complexity strongest-path beamforming scheme is employed at AP. Based on the theoretical analysis in Section 3.2, an approximate expression of average received SINR at AP can be obtained as: where the expressions of Q U,LB l and Q MUI l are given by (29) and (31), respectively. We assume that the MUI is effectively mitigated after the proposed MUI-aware sensor scheduling scheme, and as a result, the system is noise-limited (if the effect MUI is not negligible, by exploiting the analytical expression in (31), the EE optimization problem can be solved efficiently using the geometric programming based procedure developed in [40]). Thus, by substituting (29) into (33) and neglecting the low-order terms, the average received SINR can be further simplified as: where . The EE is defined as the SE divided by the AP's total energy consumption, that is: where the SE is defined as the sum of all sensors' uplink achievable rates, that is: and P c denotes the circuit power consumption of AP. According to [21,41], it is reasonable to assume that the circuit power consumption increases linearly with the number of the AP's active antennas. Thus, P c can be expressed as: P c =Ñ xÑy p ant , (37) whereÑ x andÑ y denote the numbers of active antennas along the x axis and y axis, respectively, which satisfyÑ x ≤ N x andÑ y ≤ N y . p ant denotes the power consumption introduced by each active antenna. We consider the EE optimization problem with the per sensor achievable rate constraint, which is formulated as: max where r T denotes the minimum data rate requirement of each sensor. p max denotes the maximum transmit power constraint of AP. Note that the constraint C 1 implicitly enforces the harvested energy to be greater than the processing energy. Thus, no additional constraint is needed. Due to the tremendous gain of a large-scale antenna array in a massive MIMO system, it is reasonable to assume that the uplink received SINR given by (34) is much greater than one. Thus, at a good operating point, the increase of N x , N y and p l will result in the decrease of EE since the achievable rate grows logarithmically with N x , N y and p l . In this case, the constraint C 1 will satisfy with equality at optimal N x , N y and p l . Therefore, by substituting (35) and (37) into (38) and exploiting the definition of SE, the EE optimization problem can be rewritten as: By rearranging the constraint C 1 , we can obtain: With this equality, we can rewrite the problem (39) as: With a few algebraic manipulations, we can see that the inequality constraint in (41) gives rise to a lower bound onÑ xÑy , that is: Moreover, setting the first-order derivative of the target function with respect toÑ xÑy to zero, we can obtain: Since (43) is a monotonically increasing function ofÑ xÑy , it has a single solution (denoted as N xÑy = N † ) whenÑ xÑy > 0. Therefore, the target function achieves its minimum atÑ xÑy = N † . Note that N † can be simply obtained using the method of bisection. Considering the lower bound in (42), if we discard the integer constraints onÑ x andÑ y , the optimalÑ x andÑ y satisfy: When the right-hand side of (44) is greater than N x N y , the problem is infeasible. As long as the problem is feasible,Ñ x andÑ y can be selected as arbitrary positive integers, which can make (44) satisfied.
After determining the number of active antennas, the optimal transmit power of AP can be derived directly using (40). Finally, the optimal τ can be obtained by one-dimensional search over the interval (0, 1).
From the above solutions, we can obtain two important observations on the design of system parameters.

•
First, for the above EE and widely-investigated SE optimization problems [12], we note that the optimal transmit powers derived have distinct scaling behaviors with respect to the number of the AP's (active) antennas. According to (40), the optimal transmit power that maximizes the EE reduces approximately withÑ 2 xÑ 2 y if the processing energy E p is negligible and reduces with N xÑy if the processing energy is significant. In contrast, in SE optimization, the optimal transmit power is independent of the number of AP antennas [12].

•
In EE optimization, the optimal number of active antennas, i.e.,Ñ xÑy , is a non-increasing function of transmit power constraint p max and the circuit power consumption p ant . In particular, when the circuit power consumption is dominant, (44) reduces toÑ xÑy = N LB , which becomes independent of p ant . In this case, by substituting (40) into the constraint C 2 of (38), we can see that the sum transmit power constraint of AP holds with equality. This means that, to control the total circuit power consumption, AP should use its maximum transmit power to reduce the number of active antennas.

Simulation Results and Discussion
This section evaluates the performance of the proposed beamforming scheme via MATLAB simulations. It was assumed that there are 100 sensors randomly located in the coverage area of AP. The sensors were grouped based on the MUI-aware scheduling scheme in Section 3.2. For the sake of illustration, the numbers of blocks along the x-direction and y-direction were fixed at B x = B y = 5. One sensor was picked randomly from each group, and the number of total selected sensors was L = 20 (five sensor groups were empty). The large-scale fading was modeled using the statistical propagation model for air-terrestrial transmission with low-altitude AP, which was proposed in [14]. In particular, the path loss was modeled as the summation of the free space path loss (in dB) and an extra path loss term (in dB). According to Figure 3 [14], it was assumed that the extra path loss term for each sensor varied randomly from 0 dB to 10 dB. The PAS of the channel was modeled as S (ϕ, θ) = S h (ϕ) S v (θ), where S h (ϕ) and S v (θ) denote the PASs in horizontal and vertical directions, respectively. As in [15,31], S h (ϕ) and S v (θ) were modeled using the von Mises distribution and truncated Laplacian distribution, respectively. The simulation parameters are summarized in Table 1.  We first compare the SE performance of proposed beamforming schemes with MF beamforming [10] and DFT beamforming [13] in Figures 3 and 4. Note that the MF beamforming required the instantaneous CSI and full number of RF chains (equal to N x N y ) at AP. The requirement on CSI and the number RF chains for DFT beamforming was the same as the strongest-path beamforming. Figure 4 simulates the SE performance as a function of the duration of WET phase τ. With the full number of RF chains, it was seen that the proposed statistical max-SINR beamforming achieved a similar performance with MF beamforming at the optimal τ. This makes the statistical max-SINR beamforming attractive, since it did not need the estimation of instantaneous CSI. With the limited number of RF chains, the strongest-path beamforming outperformed the DFT beamforming considerably, since it focused the signal power on the sensor's direction in a more effective way. When comparing with the MF beamforming, the strongest-path beamforming suffered from an SE loss of 11% due to the lack of instantaneous CSI.  Figure 5 simulates the SE performance for different AP transmit powers, where ∑ L l=1 p l increases from 28 dBm to 48 dBm. The duration of the WET phase was optimized with respect to SE. It is seen that the SE loss of strongest-path beamforming to the widely-used MF beamforming was quite small (about 6%) for large transmit SNR. However, the SE loss to the statistical max-SINR scheme became larger with the increase of transmit power. This was because the effect of residual MUI was dominant in the large SNR region.  Figure 6 shows the effect of sensor's processing energy E p l on SE performance, where we assumed that the processing energies for all sensors are the same and vary from 10 −11 J to 10 −8 J, as suggested in [42]. The duration of the WET phase was optimized with respect to SE. As expected, the SEs of all schemes approached zero with the increase of E p l since the energy and time resources left for the WIT phase became less. In this case, AP must increase the transmit power or equip more antennas in order to maintain the system performance. Figure 3 shows the effect of AP's altitude on the SE performance of proposed beamforming schemes. The duration of the WET phase was optimized to maximize the SE. It is seen that the SE performances of all schemes degraded with the altitude of AP since more path loss was introduced. At last, it is noted from Figures 3 and 4 that the analytical SE based on (33) gave a good approximation to the Monte Carlo result.  Figure 7 shows the EE performance of the strongest-path beamforming as a function of per antenna circuit power consumption p ant . The proposed EE optimization scheme was compared with that in [21] where the maximum number of antennas was used at AP (i.e.,Ñ x =Ñ y = 25). The figure shows that the EE performance was greatly improved by optimizing the number of the AP's active antennas, especially when the per antenna circuit power consumption p ant was large. In particular, it is seen that the performance gain increased from 1.9 times to 5.2 times when p ant varied between 10 dBm and 30 dBm.  Figure 8 simulates the effect of processing energy of sensor E p l on the proposed EE optimization scheme. It is seen that, for small E p l , the proposed optimization scheme can provide considerable gain. When E p l was greater than 10 −9 J, the EE performance under the proposed optimization scheme converged to that of [21]. This is because in this case, AP must utilize the maximum number of antennas in order to meet the achievable rate constraint and processing energy consumption of the sensor. As a result, the same solution was obtained by the two optimization schemes. Figure 9 shows the effect of AP's altitude on the proposed EE optimization scheme. Similarly, we can see that the two optimization schemes gave rise to the same EE performance for a large altitude of AP.

Conclusions
This paper studies the 3D beamforming for an airborne massive MIMO system with WET. When large numbers of RF chains are available at AP, a statistical max-SINR beamforming scheme is proposed to maximize the average received SINR at AP. A heuristic strongest-path beamforming scheme is also proposed when the number of AP RF chains is limited. Under the strongest-path beamforming scheme, the transmit power, number of active antennas of AP and duration of the WET phase are jointly optimized to maximize the system EE. The simulation results show that the statistical max-SINR beamforming scheme can achieve similar SE performance with MF beamforming. The strongest-path beamforming outperforms the DFT beamforming considerably, but suffers from a small SE loss when compared with MF beamforming. Moreover, the proposed EE optimization scheme outperforms the traditional scheme that does not consider the circuit power consumption significantly.
Funding: This work is supported by the National Natural Science Foundation of China (No. 61671472) and Jiangsu Province Natural Science Foundation (BK20160079).

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: × sin N y π ρ y − ρ y sin π ρ y − ρ y dρ x dρ y . (A3) According to [43], the function 1 N sin(Nπx) sin(πx) (this function is called aliased sinc function in [43]) converges to standard sinc function sinc(Nx) as N → ∞. Therefore, using the limit 1 a sinc x a a→0 → δ (x), we can obtain 1 N sin(Nπx) sin(πx) N→∞ → δ (x). With the above result, in the large N x , N y limit we can rewrite (A3) as 1 N x N y C NLoS l b ρ x ⊗ a ρ y = 1 N x N y β l K l + 1 S l ρ x , ρ y b ρ x ⊗ a ρ y . (A4) Note that in (A4), if ρ x / ∈ ρ min x,l , ρ max x,l or ρ y / ∈ ρ min y,l , ρ max y,l , we can simply set S l ρ x , ρ y to zero. Equation (A4) indicates that 1 √ N x N y b (ρ x ) ⊗ a ρ y is the eigenvector of C NLoS l , and the corresponding eigenvalue can be expressed as β l K l +1 S l ρ x , ρ y . For the commonly used PAS models, S l ρ x , ρ y achieves the maximum at {ρ x , ρ y } = {ρ x,l , ρ y,l } [15,36], i.e., the virtual AoDs corresponding to sensors' physical AoDs in horizontal and vertical directions. Therefore, we can conclude Combining (A2) and (A5), in the large N x , N y limit, we have u max {C l } = 1 N x N y b (ρ x,l ) ⊗ a ρ y,l , λ max (C l ) = β l K l K l + 1 N x N y + β l K l + 1 S l ρ x,l , ρ y,l . (A6)

Appendix C. Proof of Lemma 2
We first prove the first part of the lemma. Without loss of generality, we assume that only condition (a) is satisfied. From (31), the power of MUI can be rewritten as (A12) The first and second terms of right-hand side represent the MUIs contributed by the LoS and NLoS channels, respectively. By using the expression of P k in Lemma 1 and (A10), we can obtain (in the large {N x , N y } limit) 1 N x N y L ∑ k=1,k =l P k P l sin 2 (N x π (ρ x,l − ρ x,k )) N x sin 2 (π (ρ x,l − ρ x,k )) sin 2 N y π ρ y,l − ρ y,k N y sin 2 π ρ y,l − ρ y,k . (A13) Note that the function 1 N sin(Nπx) sin(πx) converges to standard sinc function sinc(Nx) as N → ∞ [43]. Thus, when condition (a) is satisfied, i.e., ρ min x,l , ρ max x,l ∩ ρ min x,k , ρ max x,k = ∅, ∀l = k, we can obtain sin 2 (N x π (ρ x,l − ρ x,k )) N x sin 2 (π (ρ x,l − ρ x,k )) = N x sin 2 (N x π (ρ x,l − ρ x,k )) where the second step follows from the inequality | sin x| ≤ 1. Moreover, the function 1 N sin(Nπx) sin (πx) achieves its maximum at point x = 0, thus, we have sin 2 N y π ρ y,l − ρ y,k N y sin 2 π ρ y,l − ρ y,k ≤ N y . When both conditions (a) and (b) are satisfied, the proof is similar with that in appendix of [44]. Thus, it is omitted.