Due to the flexibility and affordability of low-cost aerial base stations (ABSs), the coexistence of ABSs and terrestrial base stations (TBSs) equipped with a large number of antennas can be a promising solution to meet the high throughput and coverage needs of future wireless networks. In this article, we propose an analytical framework for three-dimensional (3D) heterogeneous aerial-terrestrial networks, where the TBSs are equipped with large antenna elements, and ABSs provide seamless connectivity from air to ground user equipment (GUE). Using tools from stochastic geometry, we derive exact analytical expressions for association probability, coverage probability, area spectral efficiency (ASE), and average ergodic rate (AER) while considering the maximum biased received power association scheme, path loss with line-of-sight (LoS)/non-line-of-sight (NLoS) link and Nakagami-m fading for aerial communications. The derived expressions are validated using Monte Carlo simulations, and a close match is observed between the analytical and simulated results. Simulation results demonstrate that the use of multi-antennas TBSs in conjunction with the ABSs leads to better network performance. However, the overall performance depends upon the density of the aerial and terrestrial network nodes, ABSs altitudes, and the number of antennas per TBS. Furthermore, it is observed that the dense deployment of ABSs leads to higher aerial interference. However, deploying a higher number of TBS antennas while optimally selecting other network parameters results in improved network performance.