3D Cluster-Based Ray Tracing Technique for Massive MIMO Channel Modeling

In this paper, a novel 3-dimensional (3D) approach is proposed for precise modeling of massive multiple input multiple output (M-MIMO) channels in millimeter wave (mmW) frequencies. This model is based on both deterministic and statistic computations to extract characteristics of the propagation channel. In order to increase algorithm execution speed, the physical channel is divided into two regions. The first region refers to those parts of channel which can be mapped with simple planes such as walls, ramps and etc. The second region is usually complex which is defined by considering the channel with physical clusters. These physical clusters yield multipath components (MPCs) with similar angles of arrival (AoA) and time delay. The ray tracing algorithm is utilized to find ray paths from transmitter (Tx) to receiver (Rx). Some characteristics of MPCs in each cluster are defined according to some appropriate statistical distribution. The non-stationary property of M-MIMO along the antenna array axis is considered in the algorithm. Due to the correspondence between propagation environment and scatters, the accuracy of the model is highly increased. To evaluate the proposed channel model, simulation results are compared with some measurements reported in the literature. INDEX TERMS Channel Modeling, Cluster, Massive MIMO, Stochastic Modeling.


I. INTRODUCTION
he fifth generation (5G) of wireless communication systems has attracted extensive amount of research efforts and attentions. The 5G technology can present extremely low latency, high energy efficiency, and very high data rate [1]- [5]. To achieve 5G design targets, information theory suggests three main key approaches [2,6]: (i) ultradense networks which are known as small cell technology, (ii) large quantities of new bandwidth which means utilizing higher frequencies, and (iii) high spectrum efficiency by using numerous antennas at the base station (BS) refers to massive multiple input multiple output (M-MIMO). The antenna configuration in M-MIMO systems includes hundreds or even thousands of antenna elements to increase the channel capacity, improve spectral and energy efficiency and promote reliability of the system with respect to the conventional MIMO [7]- [16]. The M-MIMO systems may confront with non-stationary properties across the array axis. This phenomenon is usually neglected in the case of conventional MIMO channels [17]. Thus, different antenna elements may face distinctive propagation characteristics. Due to the large dimension of the antenna array, a receiver may be located in a shorter distance than Rayleigh distance criteria [14]. In such a case, the far-field (FF) condition is not fulfilled, and a spherical wavefront has to be considered [15,16]. Millimeter wave technology is one of the most effective solutions to achieve huge bandwidth for 5G systems [17]- [20]. High data rates of several gigabits per second can be easily obtained in mmW communications [21]- [23].
Consequently, a mmW M-MIMO system has a great potential to improve wireless access and throughput [6,24]. Furthermore, very short wavelength of mmW frequencies helps to significantly reduce the M-MIMO array size. Additionally, the high gain of large number of antenna provided by M-MIMO systems is helpful to overcome the severe pathloss of mmW signals [25]. Since performance of mmW M-MIMO systems is greatly related to the propagation environment, channel models of mmW M-MIMO systems are very necessary in system design. Even though many models such as WINNER II [26], WINNER+ [27], COST 2100 [28], IMT-A [29] and 3D SCM [30] have been presented for conventional MIMO systems, but these channel models cannot sufficiently cover new emerging mmW M-MIMO technology for 5G systems. In recent years, many efforts have been accomplished to extract the channel behavior in mmW and M-MIMO systems.

T
The birth-death behavior of the clusters across array aperture was observed in the measurement-based channel model presented in [39], where a virtual 40×40 planar antenna array was used at Tx side. An indoor M-MIMO channel measurements at different frequencies 11, 16, 28, and 38 GHz were conducted in [40] by a large virtual uniform rectangular array. Furthermore, many geometrybased stochastic models (GBSMs) have been presented to model massive MIMO channels with new features [18], [41]- [43]. The METIS channel model is a map-based model in which the layout of the propagation environment is predefined and channel coefficients are computed based on it [43]. In [44], a 3D vehicle M-MIMO channel model has been proposed where the spherical wavefront with non-stationary conditions has been considered to statistically derive channel parameters. Another channel model has been proposed in [45], where the far-field (FF) signals are modeled by plane waves and near-field (NF) signals are modeled with spherical waves. Some cluster-based channel models are summarized in Table I. In this paper, a novel 3D hybrid channel model is proposed for wireless communication systems with the capability of simulating M-MIMO mmW channels. Here, the modeling procedure consists of two successive deterministic and stochastic modes. The combination of both modes forms the final channel model. In the deterministic mode, all surfaces, such as walls, floors, ceilings, ramps and etc. are defined with their equivalent planes. These planes are imported to the simulator in the form of rectangles along with their four corners position. Then, the Line of Sight (LoS) and reflection paths between Tx array antenna elements and Rx users of M-MIMO system are extracted, and the image theory is utilized to find the first and second order reflection paths. If a path is encountered by a surface, the transmission phenomenon is also considered in the calculations. In the stochastic mode, those parts of the channel that cannot be modeled in the deterministic mode due to their geometric complexities are defined based on the cluster concept. Accordingly, the MPCs arrive at Rx user with similar features of delay, and AoA make clusters. By means of appropriate statistical random distributions, the MPCs of the clusters are extracted and their characteristics such as phase shifts, intra-cluster delay, and their position within each cluster are modeled. The channel impulse response (CIR) is obtained from combination of two deterministic and stochastic modes. Other main characteristics of the channel such as delay spread, AoA and angle of departure (AoD) are then extracted. This proposed model not only has high accuracy, but also takes little computational time to model whole M-MIMO channel. To evaluate this model, simulation results are compared with some measurements reported in [46].
The major contributions of this paper are summarized as follows: The proposed model can be applied to a wide range of propagation channels. Since the modeling procedure is performed in two successive deterministic and stochastic modes, it can adopt to the real propagation environments, despite its simplicity of implementation. Also, the nonstationary property of the channel can be modeled along the M-MIMO antenna array in this model.
The following sections of this paper are organized as follows. Section 2 describes fundamental of the proposed channel model for non-stationary M-MIMO channels. Theoretical and mathematical details of the proposed channel model are explained in Sec.3. The scenario description is presented in Sec.4. Some simulation results and discussions are given in Sec.5. Finally, conclusions are presented in Sec.6.

II. DESCRIPTION OF THE CHANNEL MODEL
At first step, geometrical boundaries of the propagation environment as well as walls, ceiling and floor are defined and their electromagnetic properties are specified. Furthermore, Tx and Rx antenna positions and their characteristics such as radiation pattern, array elements configuration and polarization are defined. Other parameters such as frequency and transmit power are also specified. The desired channel parameters are defined at this step to be extracted at the end of the simulation procedure. At the next step, channel simulation is performed according to the predefined propagation environment and simulation strategy. To overcome the complexity of the channel, this step is divided into two deterministic and stochastic modes as mentioned before.

A. Deterministic modeling
Deterministic modeling considers those parts of the channel that can be defined as planes. If these parts are located in the FF region, the FF approximation is considered to compute the amplitude and phase shift of each path. Otherwise, the NF conditions are considered in calculations.
The 3D ray tracing algorithm is used in this stage to find propagation paths between the Tx and Rx. Two important propagation phenomena (reflection and transmission) are considered in deterministic modeling. According to Fig. 1, by utilizing image theory, the reflection points on each wall are obtained. The incident angle of each path and electromagnetic characteristics of the obstacles defined at the first step, are used to calculate Fresnel reflection and transmission coefficients according to the polarization. These reflected and transmitted paths are considered to calculate the intensity of the received signal. The whole procedures are done between each Tx element and all users. This algorithm is continued until all defined deterministic parts are checked to find existence (or non-existence) of a path between each Tx antenna and every user. Then, the algorithm runs out to switch to the stochastic modeling (if necessary).

B. STOCHASTIC MODELING
Generally, deterministic modeling of the whole channel is very complicated for M-MIMO systems. To overcome this complexity, a cluster-based stochastic modeling is proposed to complete the previous step and improve model accuracy.
The cluster-based model is shown in Fig. 2. In this model, those parts of the channel that modeled stochastically are mapped by clusters. The stochastic modeling is based on the clustering behavior of the propagation channel since the MPCs reach the Rx in the form of bundles of rays with similar properties such as delays and AoAs. The stochastic parts of the channel are defined as shown in Fig. 3a and called interacting objects (IOs). These interacting objects are fences, crowds of vehicles, vegetation and so on that cannot be modeled in the deterministic mode due to the complexity of their geometry and the lack of conditions for utilizing ray tracing algorithm. The physical boundaries of these parts are then determined, as shown in Fig. 3 b. In order to simplify the calculations, equivalent spheres are embedded in the specified area instead of IOs, which have the role of meshing in the calculations of the stochastic mode. The configuration of the equivalent spheres can have a uniform or non-uniform structure similar to that shown in Fig. 3 c. Then, the space around the Rx is divided into cluster cells as shown in Fig. 3 d. If the cluster cells have physical intersection with the equivalent spheres, then we will have an active cluster cell. Several active cluster cells are shown in Fig. 3 d. The MPCs are assigned to the active cluster cells as shown in Fig. 3 e.
The phenomena of reflection and transmission occur in the cluster and other components are neglected in the calculations due to their low amplitude level.
All MPCs experience extra phase shifts and delays in the cluster that need to be considered. The phase shifts of each MPC can be any amount in the interval of . A uniform distribution is proposed to model phase shift behavior of each MPC. The number of MPCs in each cluster is modeled by Poisson distribution with mean value of 10. An exponential distribution is used to model the intra-cluster delay time.  DOD & DOA Wrapped Gaussian or Von Mises [18], [49] The mean value of the exponential distribution is chosen according to the physical size of the clusters in which the ray strength reduced significantly and here is considered to be about 2.5 ns. These distributions are summarized in Table II.

C. COMBINATION OF TWO STAGES
As described in the previous subsections, the proposed channel model is executed in two stages. When both deterministic and statistic simulations are completed, the results should be post-processed. For this purpose, all MPCs ) 2 , 0 [ π received by the user are sorted according to their delays. The multipath delay axis is quantized into equal time delay segments called delay bins. Each bin has a time delay width equal to for to , where represents total number of equally-spaced bins, and is determined by maximum delay of MPCs.
Any MPC received within the th bin is represented by a single resolvable delay . All MPCs within delay bin of are summed up together to obtain signal strength at that bin. Then, desired characteristics of the channel are extracted based on processing of MPCs data. A general flowchart of the proposed channel model algorithm is illustrated in Fig. 4.

III. THEORETICAL DETAILS OF THE PROPOSED CHANNEL MODELING ALGORITHM
The MIMO channel is described by a matrix , where the BS is equipped with M antenna elements and serves N single-antenna users. Each entry of the matrix includes superposition of both deterministic and stochastic impulse responses. If the channel can be considered noiseless, the input-output relationship in the time domain can be expressed by: (1) where is an vector of transmitted signals, is an vector of received signal and operator (*) denotes convolution. The channel matrix is given by: where is CIR of the th time bin denotes with .

A. Mathematic theory of the deterministic algorithm
In the deterministic algorithm, the normalized unit vectors of the surfaces defined in Sec.

B. Mathematic model of stochastic algorithm
Interacting objects modeled in the stochastic mode are replaced by equivalent spheres. The configuration of the equivalent spheres is similar to the structure of the IOs. Then, active clusters are identified through the intersection detection algorithm of equivalent spheres and clusters. Although equivalent spheres can have variable centers and radii, the following limitations must also be considered: • Equivalent spheres are just defined in the area of IOs, • Depending on the distance from the Rx and the dimensions of the cluster cells, the radius of the equivalent spheres is variable, as shown in Fig. 3, • Each equivalent sphere should not occupy more than one cluster cell, • The equivalent spheres diameters should not be larger than IOs dimensions.
According to the cluster parameters shown in Fig. 5, the following condition is considered to determine the radius of each equivalent sphere :

…….where …
The range and azimuth and elevation angles of the th cluster cell are denoted by , and as shown in Fig.   5. The dimensions of the IOs are assumed larger than wavelength in the calculations. Thus, the minimum radius of the equivalent spheres is greater than half of the wavelength. The delays of each MPC are calculated based on the length of the path that each MPC travels, while the angles of arrival to the receiver depend on the position of the each MPC in the cluster. In this way, in each cluster, a number of MPCs may reach the receiver. Therefore, according to the origin of the rays generated within each cluster, the vector between the points in the corresponding cluster to the receiver is obtained. Accordingly, the arrival angles are calculated in both azimuth and elevation planes. It should be noted that the scatters are positioned based on a uniform statistical distribution in each cluster. The phase of each path is also obtained according to the path length of each MPC, in addition to the amount of phase that is added after interaction within each cluster (which is modeled using a uniform distribution). However, in relation to the amplitude of each MPC, first, the path loss is obtained according to the Friis transmission equation and then, the attenuation coefficient in each cluster is added due to the propagation phenomenon that occurs in each interaction. Since there are closed-formed mathematic relations for reflection and transmission coefficient, these coefficients are taken into account in calculating the amplitude of these MPCs. Although these two phenomena may occur in many cases, the phenomena of diffraction and scattering are also very likely. However, diffraction and scattering phenomena do not have simple and closed-form mathematical equations, and on the other hand, in most cases, the power levels of such components are very small to be neglected in the calculations.

C. Integration of deterministic and stochastic calculations
( 1)   corresponding to each bin is calculated to determine the power delay profile (PDP). Accordingly, CIR of th interval can be presented as: (7) where and are number of deterministic and stochastic MPCs with delay in the th interval , respectively. By quantizing CIR, received power from the th Tx antenna to the th user can be expressed by: where is transmit power. The mismatch and cable losses are neglected in (8).
Since AoAs of the MPCs in the respective time bin vary relative to each other, the AoA is defined individually for each MPC in the integration unit.
The root-mean-squared (RMS) delay spread can be calculated by using PDP. The RMS delay shows dispersion of delay and can be expressed as follow [46]: (9) where and are delay and received power of the th path between the th Tx antenna and the th user, respectively. Depending on the position of each antenna (assuming the phase center of each antenna which is related to the center position of that antenna), the ray tracing algorithm is applied separately to all of the antenna elements.

IV. SIMULATION ENVIRONMENT
The wave propagation simulator is developed in MATLAB software. The simulator has a main engine that consists of several sub-functions. These sub-functions calculate different parts of the channel. Figure 6 shows the general structure of the simulator program. This program includes the section of initial definitions of IOs in both deterministic and stochastic modes such as, operating frequency, transmission power, polarization of the antennas, and the section of the calculation of the arrival and departure angles of each MPC, attenuation coefficients of each MPC and its phase, and finally a section which displays the results. The elapse time of the program strongly depends on the details of the propagation environment. However, since calculating the MPCs is performed in a computational procedure and is not based on the search algorithms, modeling in this method is much faster.
All simulations are conducted for the Center Hall in JiuJiao Teaching Building, Beijing Jiaotong University, China presented in [46].  The obstacles in the Hall are assumed to be fixed without any movement. The antenna array at Tx side includes 64 elements with uniform linear configuration parallel to the ground. The carrier frequency is set to 26 GHz during our simulations. The distance between successive antenna elements is equal to half of the wavelength. Transmit power is assumed to be 1 W (0 dBW) in the simulations. All Tx and Rx antennas have isotropic radiation patterns. The simulations begin with initial parameters given in Table III

V. RESULTS AND ANALYSIS
First of all, according to Fig. 8, the walls, floor and ceiling are defined in the deterministic mode. The Tx has 64 antenna elements to cover a type of M-MIMO system. The LoS and reflection MPCs are found in this stage. The first and second order reflections are considered by applying image theory to all paths between each Tx array antenna element and Rx user. If any MPC is cut off by a surface, the transmission coefficient is also taken into account. Then, the remaining parts of the environment, such as hall seats, are defined in the stochastic mode as shown in Fig. 9. The PDP diagram between the first antenna element of the Tx array and the Rx is illustrated in Fig. 10. The measurement data extracted from the reported results in [46] is used in this figure to compare with our simulation. As illustrated by authors in [46], their measurement system includes a virtual antenna array at the Tx, while a single antenna is used at the Rx. Some more details of the measurement system are described in [46]. It can be seen that the peaks and nulls of the simulation model follow substantially the measurement results.
The RMS delay of the proposed channel model is compared in Fig. 11 with those of the measurement results reported in [46]. As it can be seen, the RMS delays fluctuate across the array antenna of M-MIMO system because the MPCs power changes along the axis of the array.       The AoAs of MPCs in the elevation and azimuth planes are shown in Fig. 13 a and Fig.13 b, respectively for the first array element. Different kinds of MPCs are shown with various signs. It can be observed in Fig. 13 a, that accumulation of MPCs in elevation plane is around 100• as it is expected. This is because of the Rx antenna is located at a higher level in the middle of the clusters. Furthermore, the MPCs within each cluster have a bit deviation in AoAs in the elevation plane. On the other hand, since the Rx antenna is located amid of the clusters, the AoAs in azimuth plane spread in a wide range of angles according to Fig. 13 b.
The AoAs on the array axis in the elevation and azimuth planes are also illustrated in Fig. 14 a and Fig. 14 b, respectively.

VI. CONCLUSION
A novel channel model for M-MIMO over mmW has been proposed in this paper for 5G networks. The proposed model is divided into two deterministic and stochastic regions. In the deterministic mode, the channel is defined for the simulator with details and by applying ray tracing algorithm, all propagation paths between Tx and Rx are determined. Then intensity and other characteristics of the received signals are obtained. But, since the M-MIMO channel is generally a very complex environment, the whole channel cannot be simulated in the deterministic mode. Thus, the second mode of the channel modeling algorithm is applied to characterize different objects located in the propagation environment. In this stage, those parts of the channel that cannot be modeled in the deterministic mode are simulated by some physical clusters through appropriate statistical distributions. This channel model has a very low computational time in comparison with full deterministic models, since it uses image theory instead of meshing the propagation environment to find the paths between Tx and Rx. In addition, it employs statistic distributions to model complex parts of the channel and reduce computational time and memory. The accuracy of the model is high, because the main structures of the propagation channel are modeled in the deterministic mode and other complex parts of the channel are mapped with corresponding clusters. This approach has been applied to an indoor environment where the field measurements have been reported. Then, simulation results have been compared with measurement results to evaluate the accuracy of the proposed algorithm. A comparison between simulation and measurement shows that total behavior of the simulated channel follows measurements results.