A Neoteric Three-Dimensional Geometry-Based Stochastic Model for Massive MIMO Fading Channels in Subway Tunnels

Abstract

Wireless mobile communication systems in subway tunnels have been widely researched these years, due to increased demand for the communication applications. As a result, an accurate model is essential to effectively evaluate the communication system performance. Thus, a neoteric three-dimensional (3D) geometry-based stochastic model (GBSM) is proposed for the massive multiple-input multiple-output (MIMO) fading channels in tunnel environment. Furthermore, the statistical properties of the channel such as space-time correlation, amplitude and phase probability density are analyzed and compared with those of the traditional two-dimensional (2D) model by numerical simulations. Finally, the ergodic capacity is investigated based on the proposed model. Numerical results show that the proposed model can describe the channel in tunnels more practically.

1. Introduction

 IN recent years, the rapid development of urban rail transit has resulted in increased demand for communication applications in subway. As a result, the wireless mobile communication systems in tunnels is extensively studied. Different from the traditional mobile communication channels on the ground, those in tunnels are confined in a limited space. Moreover, the massive multiple-input multiple-output (MIMO) technologies [1-2], where multiple antennas are eployed at both the transmitter (Tx) and receiver (Rx) [3], have been receving more and more attention in tunnel communication systems. Multiple antennas are helpful to bring huge improvements in throughput and radiated energy efficiency by focusing energy into small regions of space. Other benefits of massive MIMO include extensive use of low-power components, reduced latency and robustness against intentional jamming [4].

 In order to design an applicative communication system in tunnel environment and evaluate the system performance, an accurate channel model and a detailed knowledge of its statistical properties are indispensable. There are already many useful models adapted in different environments, such as ray-tracing models [5-6], geometry-based deterministic models (GBDMs) [7], correlation-based stochastic models (CBSMs) [8-9] and geometry-based stochastic models (GBSMs) [10-11]. Because of the research convenience for the statistical characteristics and capacity, GBSMs are most widely used. In GBSMs, the propagation characteristics of radio waves mainly depend on randomly distributed scatterers in different environment [12]. Hence, GBSMs can be adapted to different communication scenarios by modifying the random distribution of scatter clusters and the shape of the scatter area.

 There are already many available GBSMs for different communication environments. A two-dimensional (2D) model with “one-ring” or “two-ring” is proposed in [13-14] for the narrowband and wideband MIMO channels. In [15], Zheng and Nie proposed a 2D scattering model and analyzed the corresponding statistical properties. However, 2D models assume that waves travel only in the horizontal plane and neglect the vertical plane so that they cannot describe the practical three-dimensional (3D) space accurately. Therefore, in [16], the authors put forward a 3D single-ellipsoids GBSM for both outdoor and indoor environments. In [17], the 3D two-cylinder GBSM was proposed for narrowband MIMO mobile-to-mobile system. In [18], authors improved the model in [17] and developed a novel 3D one-elliptic-cylinder with two-sphere GBSM for non-isotropic MIMO vehicle-to-vehicle channels. Also, a two-sphere model [19] was provided to describe closed indoor scenarios and a two-cylinder model [12] was proposed to describe high-speed scenarios.

 The aforementioned GBSMs have not considered the train’s motion and the limited space jointly, thus, they cannot be applied directly for the particular channels in tunnels. As a result, we propose a novel theoretical 3D GBSM for the non-isotropic massive MIMO channels adapted to the tunnel environment. The proposed GBSM is combined with a horizontal cylinder model and a sphere model which is inscribed with the cylinder. Line-of-sight (LOS), single-reflection (SR) and double-reflection (DR) components are considered in this model. Based on the proposed model, some relevant comprehensive statistics are derived, i.e., the space-time correlation function (STCF), amplitude and phase probability density functions (PDFs). Then, the impact of factors on the space-time correlation are investigated and the channel capacity is analyzed.

 The remainding of this paper is organized as follows. Section 2 proposes a new 3D GBSM for massive MIMO channels in subway tunnels. In section 3, the corresponding statistical properties and channel capacity are derived. The statistics and capacity are thoroughly investigated by simulations and analyzed in section 4. Finally, the corresponding research is summarized in section 5.

2. A Neoteric 3D Massive MIMO GBSM for Channels in Subway Tunnels

 A massive MIMO wireless mobile communication system in tunnels is considered in this paper. It is supposed that the Tx and Rx are equipped with M_T and M_R number of antenna elements, respectively. Tx is stationary however, Rx is mobile. Fig. 1 illustrates our proposed GBSM with a horizontal cylinder model and a sphere model inscribed with the cylinder. The horizontal cylinder of radius R_C is used to describe the tunnel wall and the sphere of radius \(\) is used to mimic the moving trains. In Fig. 1, the antenna element at the Tx is numbered as \(p=1,2, \ldots, M_{T}\)  and the one at the Rx is numbered as \(q=1,2, \ldots, M_{R}\), respectively. The antennas at both Tx and Rx can be extended to arbitrary numbers, and the antenna arrays can be configured with any shapes. Moreover, the communication channel is depicted with the composition of the LOS, SR and DR components. The SR radio waves are incident on the cylinder and the DR ones are incident on both the sphere and cylinder. Other bounced paths are ignored because they reflect similar channel characteristics and the signal power is reduced sharply after multiple reflections [12]. In addition, there are N₁ omnidirectional effective scatterers distributed on the sphere model and N₂ scatters on the cylinder one. The n₁ th \(\left(n_{1}=1,2, \ldots, N_{1}\right)\) receive scatterer is signified by \(S^{\left(n_{1}\right)}\) and the n₂ th \(\left(n_{2}=1,2, \ldots, N_{2}\right)\) scatterer is signified by \(S^{\left(n_{2}\right)}\). The parameters in Fig. 1 are defined in Table 1. Considering the actual environment, reasonable assumptions \(D>>R_{S}\) and max, \(\max \left\{\delta_{T}, \delta_{R}\right\}<<R_{S}\) are used in the proposed GBSM.

Fig. 1. The 3D massive MIMO GBSM combining a horizontal cylinder model and a sphere one

Table 1. Definitions of symbols used in Fig. 1

 The 3D massive MIMO channel in tunnels can be described by \(\mathbf{H}(t)=\left[h_{p q}(t)\right]_{M_{T} \times M_{R}}\) matrix . The element \(h_{p q}(t)\) in the matrix denotes the received complex fading envelope of the link between the pth transmitting antenna and the qth receiving antenna at the carrier frequency of \(f_{c} \cdot h_{p q}(t)\) is a superposition of the LOS, SR and DR components, which can be expressed as [18]

\(h_{p q}(t)=h_{p q}^{\mathrm{LOS}}(t)+\sum_{i=1}^{I} h_{p q}^{\mathrm{SR}_{i}}(t)+h_{p q}^{\mathrm{DR}}(t)\)       (1)

 where I = 2 , which means there are two subcomponents for SR rays. Furthermore, SR₁ appears on the sphere and SR₂ appears on the cylinder. The complex fading envelope of each part can be respectively written as

\(h_{p q}^{\operatorname{LOS}}(t)=\sqrt{\frac{K}{K+1}} \exp \left(-\frac{\mathrm{j} 2 \pi}{\lambda} \varepsilon_{p q}\right) \cdot \exp \left[\mathrm{j} 2 \pi f_{\max } t \cos \left(\alpha_{R}^{\mathrm{LOS}}-\gamma_{R}\right) \cos \beta_{R}^{\mathrm{LOS}}\right]\)       (2)

\(\begin{aligned} h_{p q}^{\mathrm{SR}_{i}}(t) &=\sqrt{\frac{\eta_{\mathrm{SR}_{i}}}{K+1}} \lim _{N_{i} \rightarrow \infty} \sum_{n_{i}=1}^{N_{i}} \frac{1}{\sqrt{N_{i}}} \cdot \exp \left\{\mathrm{j}\left[\psi_{n_{i}}-\frac{2 \pi}{\lambda}\left(\varepsilon_{p n_{i}}+\varepsilon_{n, q}\right)\right]\right\} \\ & \cdot \exp \left\{\mathrm{j}\left[2 \pi f_{\max } t \cos \left(\alpha_{R}^{n_{i}}-\gamma_{R}\right) \cos \beta_{R}^{n_{i}}\right]\right\} \end{aligned}\)       (3)

\(\begin{aligned} h_{p q}^{\mathrm{DR}}(t) &=\sqrt{\frac{\eta_{\mathrm{DR}}}{K+1}} \lim _{N_{1}, N_{2} \rightarrow \infty} \sum_{n_{1}, n_{2}=1}^{N_{1}, N_{2}} \frac{1}{\sqrt{N_{1} N_{2}}} \cdot \exp \left\{j\left[\psi_{n_{1}, n_{2}}-\frac{2 \pi}{\lambda}\left(\varepsilon_{p n_{1}}+\varepsilon_{n_{1} n_{2}}+\varepsilon_{n_{2} q}\right)\right]\right\} \\ & \cdot \exp \left\{j\left[2 \pi f_{\max } t \cos \left(\alpha_{R}^{n_{2}}-\gamma_{R}\right) \cos \beta_{R}^{n_{2}}\right]\right\} \end{aligned}\)       (4)

where K is the Rician factor, λ is the wavelength of light, N_i is the number of effective scatterers and \(f_{\max }\) is the maximum Doppler frequency at the Rx. Here, \(\psi_{n_{i}}>\psi_{n_{1}, n_{2}}\), are independent and identically distributed (i.i.d) random variables with uniform distributions between the range \([-\pi, \pi), \quad \eta_{\mathrm{SR}_{i}}\) and \(\eta_{\mathrm{DR}}\) are power-related coefficients satisfying \(\sum_{i=1}^{I} \eta_{\mathrm{SR}_{i}}+\eta_{\mathrm{DR}}=1 . \varepsilon_{p q}, \varepsilon_{p n_{i}}, \varepsilon_{n_{i} q}\) and \(\varepsilon_{n_{1} n_{2}}\) denote \( d\left(\mathrm{Tx}_{p}, \mathrm{Rx}_{q}\right), d\left(\mathrm{Tx}_{p}, S^{\left(n_{i}\right)}\right), d\left(S^{\left(n_{i}\right)}, \mathrm{Rx}_{q}\right)\) and \(d\left(S^{\left(n_{1}\right)}, S^{\left(n_{2}\right)}\right)\) (d(X,Y) denotes the X-to-Y distance, X and Y are position points), respectively, which can be expressed as follows:

\(\begin{aligned} \varepsilon_{p q} &=\left[\left(D-\frac{\delta_{T}}{2} \cos \varphi_{T} \cos \theta_{T}+\frac{\delta_{R}}{2} \cos \varphi_{R} \cos \theta_{R}\right)^{2}\right.\\ &\left.+\left(\frac{\delta_{T}}{2} \cos \varphi_{T} \sin \theta_{T}-\frac{\delta_{R}}{2} \cos \varphi_{R} \sin \theta_{R}\right)^{2}+\left(\frac{\delta_{T}}{2} \sin \varphi_{T}-\frac{\delta_{R}}{2} \sin \varphi_{R}\right)^{2}\right]^{1 / 2} \end{aligned}\)       (5)

\(\begin{aligned} \varepsilon_{p n_{1}} &=\left[\left(\frac{\delta_{T}}{2} \cos \varphi_{T} \sin \theta_{T}-R_{S} \cos \beta_{R}^{n_{1}} \sin \alpha_{R}^{n_{1}}\right)^{2}\right.\\ &\left.+\left(D-\frac{\delta_{T}}{2} \cos \varphi_{T} \cos \theta_{T}+R_{S} \cos \beta_{R}^{n_{1}} \cos \alpha_{R}^{n_{1}}\right)^{2}+\left(\frac{\delta_{T}}{2} \sin \varphi_{T}-R_{S} \sin \beta_{R}^{n_{1}}\right)^{2}\right]^{1 / 2} \end{aligned}\)       (6)

\(\begin{aligned} \varepsilon_{n_{1} q} &=\left[\left(\frac{\delta_{R}}{2} \cos \varphi_{R} \cos \theta_{R}-R_{S} \cos \beta_{R}^{n_{1}} \cos \alpha_{R}^{n_{1}}\right)^{2}\right.\\ &\left.+\left(\frac{\delta_{R}}{2} \cos \varphi_{R} \sin \theta_{R}-R_{S} \cos \beta_{R}^{n_{1}} \sin \alpha_{R}^{n_{1}}\right)^{2}+\left(\frac{\delta_{R}}{2} \sin \varphi_{R}-R_{S} \sin \beta_{R}^{n_{1}}\right)^{2}\right]^{1 / 2} \end{aligned}\)       (7)

\(\begin{aligned} \varepsilon_{p n_{2}} &=\left[\left(\xi_{T{n_{2}}} \cos \beta_{T}^{n_{2}} \cos \alpha_{T}^{n_{2}}-\frac{\delta_{T}}{2} \cos \varphi_{T} \cos \theta_{T}\right)^{2}\right.\\ &\left.+\left(\frac{\delta_{T}}{2} \cos \varphi_{T} \sin \theta_{T}-\xi_{T n_{2}} \cos \beta_{T}^{n_{2}} \sin \alpha_{T}^{n_{2}}\right)^{2}+\left(\frac{\delta_{T}}{2} \sin \varphi_{T}+\xi_{T{n_{2}}} \sin \beta_{T}^{n_{2}}\right)^{2}\right]^{1 / 2} \end{aligned}\)       (8)

\(\begin{aligned} \varepsilon_{n_{2} q} &=\left[\left(\frac{\delta_{R}}{2} \cos \varphi_{R} \cos \theta_{R}-\xi_{R n_{2}} \cos \beta_{R}^{n_{2}} \cos \alpha_{R}^{n_{2}}\right)^{2}\right.\\ &\left.+\left(\frac{\delta_{R}}{2} \cos \varphi_{R} \sin \theta_{R}-\xi_{R n_{2}} \cos \beta_{R}^{n_{2}} \sin \alpha_{R}^{n_{2}}\right)^{2}+\left(\frac{\delta_{R}}{2} \sin \varphi_{R}+\xi_{R n_{2}} \sin \beta_{R}^{n_{2}}\right)^{2}\right]^{1 / 2} \end{aligned}\)       (9)

\(\begin{aligned} \varepsilon_{n_{1} n_{2}} &=\left[\left(R_{S} \sin \beta_{R}^{n_{1}}+\xi_{R n_{2}} \sin \beta_{R}^{n_{2}}\right)^{2}+\left(R_{S} \cos \beta_{R}^{n_{1}} \cos \alpha_{R}^{n_{1}}-\xi_{R n_{2}} \cos \beta_{R}^{n_{2}} \cos \alpha_{R}^{n_{2}}\right)^{2}\right.\\ &\left.+\left(R_{S} \cos \beta_{R}^{n_{1}} \sin \alpha_{R}^{n_{1}}-\xi_{R n_{2}} \cos \beta_{R}^{n_{2}} \sin \alpha_{R}^{n_{2}}\right)^{2}\right]^{1 / 2} \end{aligned}\)       (10)

where \(\xi_{T n_{2}}=R_{C} / \sin \left[\arccos \left(\cos \beta_{T}^{n_{2}} \cos \alpha_{T}^{n_{2}}\right)\right]\) and \(\xi_{R n_{2}}=R_{C} / \sin \left[\arccos \left(\cos \beta_{R}^{n_{2}} \cos \left(\pi-\alpha_{R}^{n_{2}}\right)\right)\right]\).

 Note that the determined scatterer distribution needs to be applied in this particular model for the subsequent study of relevant statistical properties. In this paper, to consider the effect of the azimuth and elevation angles jointly, von Mises-Fisher (VMF) distribution is used to characterize the effective scatterer distribution, whose PDF is defined as [20]

\(f(\alpha, \beta)=\frac{k \cos \beta}{4 \pi \sinh k} \cdot \exp \left[k\left(\cos \beta \cos \beta_{0} \cos \left(\alpha-\alpha_{0}\right)+\sin \beta \sin \beta_{0}\right)\right]\)       (11)

where \(\alpha \in[-\pi, \pi), \beta \in[-\pi / 2, \pi / 2), \alpha_{0} \in[-\pi, \pi)\) and \(\beta_{0} \in[-\pi / 2, \pi / 2)\) represent the average values of azimuth and elevation angles, respectively, and k is an environmental coefficient which controls the scatterer spread around the mean value.

3. Statistical Properties and Channel Capacity of the 3D Massive MIMO GBSM

 The massive MIMO system is capable of making full use of space resources and can achieve the spatial parallel transmission of radio waves. Furthermore, the link reliability and channel capacity can be improved significantly by using space-time joint processing method [21]. In order to evaluate the performance of the 3D GBSM in tunnels more effectively and associate the effective scatterers in the propagation environment with the multiple antennas’ performance, the statistical properties, such as STCF and PDFs, and the capacity need to be analyzed.

3.1 STCF

 Both the spatial and temporal correlations exist in any fading channels. The spatial correlation describes the cross-correlation properties between different antenna elements while the temporal one depicts the auto-correlation properties of a single receiving antenna element at different time instances. The STCF can be used to describe the impact of both spatial and temporal correlations on wireless communication systems.

 Under the wide sense stationary uncorrelated scattering condition [22], the normalized STCF between any two sub-channels can be derived from the complex fading envelope showed in (1) and written as [23]

\(\begin{aligned} \rho_{p q, p' q'}\left(\delta_{T}, \delta_{R}, \tau\right) &=\frac{\mathbf{E}\left[h_{p q}(t) h_{p' q'}^{*}(t+\tau)\right]}{\sqrt{\mathbf{E}\left[\left|h_{p q}(t)\right|^{2}\right] \cdot \mathbf{E}\left[\left|h_{p^{\prime} q'}(t)\right|^{2}\right]}} \\ &=\rho_{p q, p' q'}^{\mathrm{LOS}}\left(\delta_{T}, \delta_{R}, \tau\right)+\sum_{i=1}^{I} \rho_{p q, p' q'}^{\mathrm{SR}_{i}}\left(\delta_{T}, \delta_{R}, \tau\right)+\rho_{p q, p' q'}^{\mathrm{DR}}\left(\delta_{T}, \delta_{R}, \tau\right) \end{aligned}\)       (12)

where E[⋅] denotes the statistical expectation operator and (⋅)^∗ is the complex conjugate operation, \(p, p^{\prime} \in\left(1,2, \ldots, M_{T}\right)\)) and \(q, q^{\prime} \in\left(1,2, \ldots, M_{R}\right)\).

 Substituting (2)-(4) into (12) and using the accordant VMF scattering distribution at the same time, the STCF of LOS, SR and DR components is derived next.

 1) The STCF of LOS component

 The STCF of LOS component is expressed as

\(\rho_{p q, p' q'}^{\mathrm{LOS}}\left(\delta_{T}, \delta_{R}, \tau\right)=K \cdot \exp \left(\frac{\mathrm{j} 2 \pi}{\lambda} \cdot D_{\mathrm{LOS}}\right) \cdot \exp \left[\mathrm{j} 2 \pi \tau f_{\max } \cos \left(\alpha_{R}^{\mathrm{LOS}}-\gamma_{R}\right) \cos \beta_{R}^{\mathrm{LOS}}\right]\)       (13)

where \(D_{\mathrm{LOS}}=2 D \cos \varphi_{R} \cos \theta_{R}\). Furthermore, considering D >> R_S , we have the following approximations \(\alpha_{R}^{\mathrm{LOS}} \approx \pi\) and \(\alpha_{T}^{\mathrm{LOS}} \approx \beta_{T}^{\mathrm{LOS}} \approx \beta_{R}^{\mathrm{LOS}} \approx 0\). Therefore, Eq. (13) can be simplified as

\(\rho_{p q, p' q'}^{\mathrm{LOS}}\left(\delta_{T}, \delta_{R}, \tau\right)=K \cdot \exp \left(\frac{\mathrm{j} 2 \pi}{\lambda} D_{\mathrm{LOS}}\right) \cdot \exp \left(-\mathrm{j} 2 \pi \tau f_{\max } \cos \gamma_{R}\right)\)       (14)

 2) The STCF of SR component

 The STCF of SR component is expressed as

\(\begin{aligned} \rho_{p q, p' q'}^{\mathrm{SR}_{i}}\left(\delta_{T}, \delta_{R}, \tau\right) &=\eta_{\mathrm{SR}_{i}} \int_{-\pi}^{\pi} \int_{-\pi}^{\pi}\left\{\frac{k_{i} \cos \beta_{R_{i}}^{n_{i}}}{4 \pi \sinh k_{i}} \cdot \exp \left(\frac{\mathrm{j} 2 \pi D_{\mathrm{SR}_{i}}}{\lambda}\right)\right.\\ & \cdot \exp \left(\mathrm{j} 2 \pi \tau f_{\max } \cos \left(\alpha_{R}^{n_{i}}-\gamma_{R}\right) \cos \beta_{R}^{n_{i}}\right) \\ &\left.\cdot \exp \left[k_{i}\left(\cos \beta_{R}^{n_{i}} \cos \beta_{R 0}^{n_{i}} \cos \left(\alpha_{R}^{n_{i}}-\alpha_{R 0}^{n_{i}}\right)+\sin \beta_{R}^{n_{i}} \sin \beta_{R 0}^{n_{i}}\right)\right]\right\} d \alpha_{R}^{n_{i}} d \beta_{R}^{n_{i}} \end{aligned}\)       (15)

where \(D_{\mathrm{SR}_{1}}=\delta_{T}\left[R_{S} \sin \beta_{R}^{n_{1}} \sin \varphi_{T}+Q_{n_{1}} \cos \varphi_{T} \cos \left(\theta_{T}-\alpha_{T}^{n_{1}}\right)\right] / \xi_{T n_{1}}+\delta_{R}\left[\sin \beta_{R}^{n_{1}} \sin \varphi_{R}\right.\)\(\left.+\cos \beta_{R}^{n_{1}} \cos \varphi_{R} \cos \left(\theta_{R}-\alpha_{R}^{n_{1}}\right)\right]\), \(D_{\mathrm{SR}_{2}}=\delta_{R}\left[\sin \beta_{R}^{n_{2}} \sin \varphi_{R}+\cos \beta_{R}^{n_{2}} \cos \varphi_{R} \cos \left(\theta_{R}-\alpha_{R}^{n_{2}}\right)\right]\)\(+\delta_{T}\left[\xi_{R n_{2}} \sin \beta_{R}^{n_{2}} \sin \varphi_{T}+Q_{n_{2}} \cos \varphi_{T} \cos \left(\theta_{T}-\alpha_{T}^{n_{2}}\right)\right] / \xi_{T_{n_{2}}}\), \(, Q_{n_{1}}=D+R_{S} \cos \beta_{R}^{n_{1}} \cos \alpha_{R}^{n_{1}}\), \(Q_{n_{2}}=\left(\xi_{T n_{2}}^{2}-\xi_{R n_{2}}^{2} \sin ^{2} \beta_{R}^{n_{2}}\right)^{1 / 2}\), \(\xi_{T n_{1}}=\left[\left(D+R_{S} \cos \beta_{R}^{n_{1}} \cos \alpha_{R}^{n_{1}}\right)^{2}+\left(R_{S} \sin \beta_{R}^{n_{1}}\right)^{2}\right]^{1 / 2}\), \(\xi_{T n_{2}}=R_{C} / \sin \left[\arccos \left(\cos \beta_{T}^{n_{2}} \cos \alpha_{T}^{n_{2}}\right)\right]\), \(\xi_{R n_{2}}=R_{C} / \sin \left[\arccos \left(\cos \beta_{R}^{n_{2}} \cos \left(\pi-\alpha_{R}^{n_{2}}\right)\right)\right]\). Due to D>>R_S and considering the relationship between the angles of the Tx and Rx, we have the following approximations such as \(\alpha_{T}^{n_{1}} \approx \beta_{T}^{n_{1}} \approx 0, \quad \xi_{T n_{1}} \approx Q_{n_{1}} \approx D\) and \(\xi_{T n_{2}} \approx \xi_{R n_{2}}\). Thus, \(D_{\mathrm{SR}_{1}} \approx \delta_{R}\left[\sin \beta_{R}^{n_{1}} \sin \varphi_{R}+\cos \beta_{R}^{n_{1}} \cos \varphi_{R} \cos \left(\theta_{R}-\alpha_{R}^{n_{1}}\right)\right]+\delta_{T}\left[R_{S} \sin \beta_{R}^{n_{1}} \sin \varphi_{T} / D+\cos \varphi_{T} \cos \theta_{T}\right]\) , \(D_{\mathrm{SR}_{2}} \approx \delta_{R}\left[\cos \beta_{R}^{n_{2}} \cos \varphi_{R} \cos \left(\theta_{R}-\alpha_{R}^{n_{2}}\right)+\sin \beta_{R}^{n_{2}} \sin \varphi_{R}\right]+\delta_{T}\left[\sin \beta_{R}^{n_{2}} \sin \varphi_{T}\right.\) \(\left.+\cos \beta_{R}^{n_{2}} \cos \varphi_{T} \cos \left(\theta_{T}-\alpha_{T}^{n_{2}}\right)\right]\).

 3) The STCF of DR component

 The STCF of DR component is expressed as

\(\begin{aligned} \rho_{p q, p' q'}^{\mathrm{DR}},\left(\delta_{T}, \delta_{R}, \tau\right) &=\eta_{\mathrm{DR}} \int_{-\pi}^{\pi} \int_{-\pi}^{\pi} \int_{-\pi}^{\pi} \int_{-\pi}^{\pi}\left\{\frac{k_{1} k_{2} \cos \beta_{R}^{n_{1}} \cos \beta_{T}^{n_{2}}}{16 \pi^{2} \sinh k_{1} \sinh k_{2}} \cdot \exp \left(\frac{\mathrm{j} 2 \pi}{\lambda} D_{\mathrm{DR}}\right)\right.\\ & \cdot \exp \left[k_{2}\left(\cos \beta_{T}^{n_{2}} \cos \beta_{T 0}^{n_{2}} \cos \left(\alpha_{T}^{n_{2}}-\alpha_{T 0}^{n_{2}}\right)+\sin \beta_{T}^{n_{2}} \sin \beta_{T 0}^{n_{2}}\right)\right] \\ & \cdot \exp \left[k_{1}\left(\cos \beta_{R}^{n_{1}} \cos \beta_{R 0}^{n_{1}} \cos \left(\alpha_{R}^{n_{1}}-\alpha_{R 0}^{n_{1}}\right)+\sin \beta_{R}^{n_{1}} \sin \beta_{R 0}^{n_{1}}\right)\right] \\ &\left.\cdot \exp \left[j 2 \pi \tau f_{\max } \cos \left(\alpha_{R}^{n_{2}}-\gamma_{R}\right) \cos \beta_{R}^{n_{2}}\right]\right\} d \alpha_{T}^{n_{2}} d \beta_{T}^{n_{2}} d \alpha_{R}^{n_{1}} d \beta_{R}^{n_{1}} \end{aligned}\)       (16)

where\(D_{\mathrm{DR}}=\delta_{T}\left[\xi_{R n_{2}} \sin \beta_{R}^{n_{2}} \sin \varphi_{T}+Q_{n_{2}} \cos \varphi_{T} \cos \left(\theta_{T}-\alpha_{T}^{n_{2}}\right)\right] / \xi_{T n_{2}}+\delta_{R}\left[\sin \beta_{R}^{n_{1}} \sin \varphi_{R}\right.\) \(\left.+\cos \beta_{R}^{n_{1}} \cos \varphi_{R} \cos \left(\theta_{R}-\alpha_{R}^{n_{1}}\right)\right]\).

3.2 Amplitude and Phase PDFs

 The amplitude and phase of the complex fading envelope \(h_{p q}(t)\)\(\) are \(z(t)=\left|h_{p q}(t)\right|\) and \(\theta(t)=\arg \left\{h_{p q}(t)\right\}\), respectively. Hence, according to the similar computation procedure offered in [24], the amplitude PDF of the 3D GBSM can be derived as

\(\begin{array}{c} g_{z}(z)=4 \pi^{2} z \int_{0}^{\infty}\left[\prod_{n_{1}=1}^{N_{1}} J_{0}\left(2 \pi\left|Q_{\mathrm{SR}_{1}}\right| x\right) \cdot \prod_{n_{2}=1}^{N_{2}} J_{0}\left(2 \pi\left|Q_{\mathrm{SR}_{2}}\right| x\right) \cdot \prod_{n_{1}, n_{2}=1}^{N_{1}, N_{2}} J_{0}\left(2 \pi\left|Q_{\mathrm{DR}}\right| x\right)\right] \\ \quad J_{0}(2 \pi z x) \cdot J_{0}\left(2 \pi K_{0} x\right) x \mathrm{d} x \end{array}\)       (17)

and the phase PDF can be similarly expressed as

\(\begin{aligned} g_{\theta}(\theta)=& 2 \pi \int_{0}^{\infty} \int_{0}^{\infty}\left[\prod_{n_{1}=1}^{N_{1}} J_{0}\left(2 \pi\left|Q_{\mathrm{SR}_{1}}\right| x\right) \cdot \prod_{n_{2}=1}^{N_{2}} J_{0}\left(2 \pi\left|Q_{\mathrm{SR}_{2}}\right| x\right) \cdot \prod_{n_{1}, n_{2}=1}^{N_{1}, N_{2}} J_{0}\left(2 \pi\left|Q_{\mathrm{DR}}\right| x\right)\right] \\ & \cdot J_{0}(2 \pi x \sqrt{z^{2}+K_{0}^{2}-2 z K_{0} \cos \left(\theta-\theta_{r e}\right)}) x z \mathrm{d} x \mathrm{d} z \end{aligned}\)       (18)

where \(Q_{\mathrm{SR}_{i}}=\sqrt{\frac{\eta_{\mathrm{SR}_{i}}}{N_{i}(K+1)}}\), \(Q_{\mathrm{DR}}=\sqrt{\frac{\eta_{\mathrm{DR}}}{N_{1} N_{2}(K+1)}} \quad( i=1,2)\),\(K_{0}=\sqrt{\frac{K}{K+1}}\), \(\theta_{r e}=\arg \left\{h_{p q}^{\mathrm{LOS}}(t)\right\}\) and \(J_{0}(\cdot)\) represents the zeroth-order Bessel function of the first kind.

3.3 Channel Capacity

In this section, the MIMO channel capacity is derived. Channel capacity refers to the maximum average information transmitting rate of a channel. Furthermore, explicit MIMO system structure and theory are proposed in literature [25] and the authors proved mathematically that channel capacity can be increased by increasing the transmitting and receiving antennas without augmenting the bandwidth. Under an average transmit power constraint, the instantaneous channel capacity of a stochastic massive MIMO channel can be defined as [26]

\(C(t)=\log _{2}\left[\operatorname{det}\left(\mathbf{I}_{\min }+\frac{\rho_{S N R}}{M_{T}} \mathbf{H}(t) \mathbf{H}^{\mathrm{H}}(t)\right)\right]\)       (19)

where ρ_SNR is the average signal-to-noise ratio(SNR) at the Tx, I_min represents the \(M_{R} \times M_{R}\) identity matrix, \((\cdot)^{\mathrm{H}}\) denotes the transpose conjugate operation and det(⋅) designates the matrix determinant. Here, it’s assumed that \(M_{R} \leq M_{T}\) which is the condition that space-division multiplexing must meet [27].

 In general, channel capacity has two significant definitions which are ergodic capacity and outage capacity. Furthermore, the ergodic capacity is the mathematic expectation of the instantaneous capacity over time and is defined as[27]

\(C_{m e a n}=\mathbf{E}[C(t)]=\mathbf{E}\left\{\log _{2}\left[\operatorname{det}\left(\mathbf{I}_{\min }+\frac{\rho_{S N R}}{M_{T}} \mathbf{H}(t) \mathbf{H}^{\mathrm{H}}(t)\right)\right]\right\}\)       (20)

and the outage capacity C_out means the most suitable capacity under the condition that message is correctly received. The C_out is related to an outage probability P_out which satisfies

\(\begin{aligned} P_{o u t} &=p\left(\mathbf{H}: C(t)<C_{o u t}\right) \\ &=p\left(\mathbf{H}: \log _{2}\left[\operatorname{det}\left(\mathbf{I}_{\min }+\frac{\rho_{S N R}}{M_{T}} \mathbf{H}(t) \mathbf{H}^{\mathrm{H}}(t)\right]\right]<C_{o u t}\right) \end{aligned}\)       (21)

where p(⋅) is a probability function.

4. Simulation Results and Analysis

 The statistical properties of the proposed 3D GBSM are investigated in this section. Referring to the literature [28-29] and the 3rd Generation Partnership Project (3GPP) protocol [30], the simulation parameters chosen are: \(f_{c}=2 \mathrm{GHz}, D=500 \mathrm{m}, R_{C}=R_{S}=4 \mathrm{m}, f_{\mathrm{max}}=500 \mathrm{Hz}, \gamma_{R}=0^{\circ}\), \(\theta_{T}=\theta_{R}=45^{\circ}, \eta_{\mathrm{SR}_{\mathrm{I}}}=0.125, \eta_{\mathrm{SR}_{2}}=0.325, \eta_{\mathrm{DR}}=0.55\). The environmental coefficient k_i is related to the effective scatterer distribution. Specifically, for the denser environmental scatterers, the smaller value of k_i is used. We set k₁=4 , k₂=15 , M_T=128 and M_R=64 . In addition, when all elevation angles are set as 0 , the 3D GBSM will be translated into 2D model which can be taken as a contrast. The performances are compared with the 3D model, 2D model and the geometric elliptical scattering model proposed in [31].

4.1 Temporal Auto-Correlation Function (ACF)

 The expression of the temporal ACF can be derived from STCF by setting \(\delta_{T}=\delta_{R}=0\) and is expressed as

\(\rho_{p q, p' q'}(\tau)=\rho_{p q, p' q'} (0,0, \tau)\)       (22)

 Fig. 2 shows the absolute values of the normalized temporal ACFs for the proposed 3D model, the 3D reference model and the 2D model. It is clear that the absolute values of ACFs decrease as the time separation increases whether in the 3D or 2D model which implies that auto-correlation values drop. Moreover, the ACF value of the 3D reference model is always lower than that of the 2D model but higher than that of the 3D proposed model.

Fig. 2. The absolute values of the temporal ACFs for the 3D and 2D models (\(\delta_{T}=\delta_{R}=0\)).

4.2 Spatial Cross-Correlation Function (CCF)

 The CCF can be used to describe the mutual influence between different antenna elements under the same time. In order to analyze the CCF, τ is set as 0 and δ_T is set as 0.5λ , thus, the spatial CCF can be expressed as

\(\rho_{p q, 'p q'}\left(\delta_{R}\right)=\rho_{p q, p' q'}\left(0.5 \lambda, \delta_{R}, 0\right)\)       (23)

 Fig. 3 presents the absolute values of the normalized spatial CCFs for the 3D model ( \(\varphi_{T}=45^{\circ}, \varphi_{R}=45^{\circ}\) and \(\varphi_{T}=45^{\circ}, \varphi_{R}=60^{\circ}\)), the 3D reference model and the 2D model (\(\varphi_{T}=\varphi_{R}=0^{\circ}\)). Similarly, it is clear that as the normalized spacing between antenna elements at the Rx increase, the absolute values of CCFs decrease, which signifies that the spatial correlation gradually diminishes. Furthermore, the CCF value of the 2D model is always higher than that of the 3D proposed model, which means the 2D model overestimates the spatial correlation. Moreover, Compared with the 3D reference model, the model in this paper has lower CCF. In addition, the elevation angle of the Rx antenna array ϕ_R has significant impact on CCFs. Specifically, when \(0 \leq \varphi_{R} \leq 90^{\circ}\), the bigger ϕ_R , the lower CCF values.

Fig. 3. The absolute values of the spatial CCFs for the 3D and 2D models (\(\tau=0, \delta_{T} / \lambda=0.5\)).

4.3 Amplitude and Phase PDFs

 Figs. 4 and 5 show that the Rician factor K has a significant influence on the amplitude and phase PDFs. Specifically, the curve of the amplitude PDF shifts to the right and the peak value decreases as K increases. As to the phase PDF, it is clear that increasing the K will improve the peak value however, the mean value remain constant.

Fig. 4. The amplitude PDFs with K=1, K=3.5 and K=9.

Fig. 5. The phase PDFs with K=1, K=3.5 and K=9.

4.4 Channel Capacity Analysis

 The factors effecting the channel capacity such as SNR and space-time correlation are investigated in this section. Fig. 6 shows the ergodic capacity against ρ_SNR with \(\delta_{T} / \lambda=\delta_{R} / \lambda=5\) for both 3D and 2D models. As expected, the ergodic capacity increases as ρ_SNR increases however, the growth rate gradually decreases. Furthermore, similar to the previous temporal and spatial correlations, the 3D model always has a higher ergodic capacity than the 2D model. Fig. 7 presents the ergodic capacity as a function of the normalized antenna element spacing at the Tx and Rx with ρ_SNR = 10 dB for the 3D model. It is observed that when \(\delta_{T} / \lambda<3\), increasing the normalized antenna element spacing \(\delta_{T} / \lambda\) increases the ergodic capacity. However, it is clear that increasing \(\delta_{T} / \lambda\) has a negligible effect on the ergodic capacity when \(\delta_{T} / \lambda>3\). In addition, the ergodic capacity can be improved by incrasing \(\delta_{R} / \lambda\) but cannot be increased indefinitely.

Fig. 6. The ergodic capacity against ρ_SNR for the 3D and 2D models.

Fig. 7. The ergodic capacity as a function of the normalized antenna element spacing at the Tx and Rx for the 3D model

5. Conclusion

In this paper, a neoteric 3D “cylinder-sphere” theoretical geometrical model for the non-isotropic scattering massive MIMO channels in tunnels is proposed. Furthermore, the matrix of complex fading envelopes for the massive MIMO channel is derived. Then, the STCF, amplitude and phase PDFs are investigated and the derived channel capacity is analyzed by numerical simulations. The statistics of the 3D and 2D models are also verified by simulations. It is obvious that the 3D model has a lower space-time correlation compared with the 2D model. Besides, the time separation and the normalized antenna element spacing at the Rx both have negative effect on the STCF and the Rician factor K has a significant impact on the amplitude and phase PDFs. Furthermore, numerical simulations have shown that the channel capacity increases as the SNR and normalized antenna element spacing increase, respectively. In a nutshell, the proposed 3D massive MIMO GBSM can characterize the channel environment in tunnels more practically than the traditional 2D models. The presented research work can be taken as a reference to establish the future communication systems operating in subway tunnels.