Analytical Nonstationary 3D MIMO Channel Model for Vehicle-to-Vehicle Communication on Slope

Vehicle-to-vehicle communication plays a strong role in modern wireless communication systems, appropriate channel models are of great importance in future research, and propagation environment with slope is one special kind. In this study, a novel three-dimensional nonstationary multiple-input multiple-output channel model for the sub-6GHz band is proposed. /is model is a regular-shaped multicluster geometry-based analytical model, and it combines the line-of-sight component and multicluster scattering rays as the nonline-of-sight components. Each cluster of scatterers represents the influence of different moving vehicles on or near a slope, and scatterers are, respectively, distributed within two spheres around the transmitter and the receiver. In this model, it is considered that the azimuth and elevation angles of departure and arrival are jointly distributed and conform to the von Mises–Fisher distribution, which can easily control the range and concentration of the scatterers within spheres to mimic the real-world situation well. Moreover, the impulse response and the autocorrelation function of the corresponding channel is derived and proposed; then, the Doppler power spectrum density of the channel is simulated and analyzed. In addition, the nonstationary characteristics of the presented channel model are observed through simulations. Finally, the simulation results are compared with measurement data in order to validate the utility of the proposed model.


Introduction
In recent years, the applications of vehicle-to-vehicle (V2V) communications have become more extensive and valuable, such as vehicular networking [1][2][3], the cooperative vehicular system [4], and intelligent transportation [5]. Besides, the accuracy of performance evaluation for V2V communications is to some extent determined by precise modeling of the channel. Reference [6] analyzed vehicular communications from the physical layer. With the mature multipleinput multiple-output (MIMO) technology [7,8], the performance of V2V communications has been greatly improved. Compared to other communication scenarios, V2V communications have their unique characteristics, i.e., both the transmitter (Tx) and the receiver (Rx) are in high-speed movements, the shapes and heights of different vehicles are similar, and the heights of vehicles are generally between one and two meters, so vehicles usually have antennas with low elevation angles. us, the obstruction of the slope will affect the performance of V2V communication seriously. With these special features, the channel modeling of V2V communications is essential and necessary, especially the propagation model on slope terrain. e path loss for V2V scenario with slope has been studied in several works. References [9,10] proposed path loss prediction models for slope terrain in small-sized cells. However, research on the modeling of scatterers still needs to be finished. e modeling of V2V channels can be of great help in testing the performance and properties of V2V channels. Reference [11] presented propagation models for urban environment. Besides, several precise channel models for other scenarios have been proposed by other researchers. According to this literature, channel models can generally be classified into deterministic models and stochastic models. Stochastic models are further categorized as geometry-based stochastic models (GBSMs) [12] and nongeometry-based stochastic models (NGBSMs) [13]. GBSMs assume that scatterers are distributed on a certain geometry and use methods similar to ray-tracing to represent scattering rays. In further classification, GBSMs can be categorized as regular-shaped GBSMs (RS-GBSMs) [14] and irregularshaped GBSMs (IS-GBSMs) [15]. RS-GBSMs assume that these geometries are regular shapes, e.g., sphere, ring, or cylinder. In other ways, GBSMs can be divided into twodimensional (2D) GBSMs [16] and three-dimensional (3D) GBSMs [17] according to the dimension of models. Considering the pattern of antenna arrays, the models can be categorized as single-input single-output (SISO) models [18] and MIMO models [19]. Additionally, in [20], a 3D model adopting massive-MIMO technology was proposed.
Comparisons with measurements [21][22][23] showed that GBSMs could be applied to V2V channel modeling suitably. V2V channels were modeled in 2D due to the very narrow elevation angle of the vehicle antennas at the early time. In [24], the authors proposed a nonstationary 2D RS-GBSM consisting of multiple confocal ellipses for high-mobility intelligent transportation systems. In [25], a 2D RS-GBSM with two-ring for a nonisotropic scattering environment was proposed. Furthermore, 2D models were extended to 3D models because they could better simulate the realistic propagation environment, e.g., the 2D two-circular RS-GBSM [25] was updated to the 3D two-cylinder RS-GBSM in [26], and the 2D two-circular and one-ellipse RS-GBSM [27] was extended to the 3D two-sphere and one-ellipse-cylinder RS-GBSM in [28]. Including the above listed, 3D GBSMs can be classified according to the geometries of scattering, such as spheres, cylinders, one/two semiellipsoid [29][30][31], and one/two hollow semiellipsoid around the Tx/Rx [32,33]. In addition, GBSMs can be divided into stationary GBSMs and nonstationary GBSMs [34][35][36][37][38][39][40][41][42]. e nonstationary GBSMs are closer to reality due to the time-varying parameters of the V2V channel. According to the abovementioned literature, it can be noticed that in the RS-GBSM based channel modeling, the modeling of the scatterer distribution is crucial. In RS-GBSM-based models, isotropic and nonisotropic distributions of scatterers are considered, and the latter is more reasonable. In [27,35,36,43], von Mises distributions were adopted to model the distributions of the scatterers in 2D modeling as nonisotropic models. In [28,37,[44][45][46][47][48], 3D nonisotropic distributions of scatterers were modeled with von Mises-Fisher (VMF) distributions. By now, the VMF distribution is perhaps the most suitable distribution for 3D nonisotropic scattering modeling.
In the aforementioned literature, although nonisotropic 3D scatterers have been modeled, scatterers were always assumed to be concentrated in the direction of the line-ofsight (LoS) path. Further discussions on multidirection scatterers modeling should be conducted. Besides, limited specific models are for vehicles near or on a slope as well as overtaking cases. e motivation of this research is to present a flexible 3D nonstationary V2V channel model with multidirection distributing scatterers on or near slope terrain. Based on the two-sphere 3D RS-GBSM V2V nonstationary models in [28,37,48,49], we abandoned the additional elliptical cylinder and introduced multidirection scattering clusters to model the multidirectional isotropic scatterers. More distinctly, in the model we proposed, rays are scattered by other vehicles in several directions in the form of cluster around Tx and Rx. is is the major difference between our model and [28,49,50] and basic theories found in [51]. For the modeling of nonisotropic scattering, the elevation and azimuth angles of arrival (AoAs) and the elevation and azimuth angles of departure (AoDs) are jointly distributed and subjected to VMF distribution, as in [45,47]. Finally, we compare the simulation results with the measurement data.
In summary, the main contributions of this study are as follows: e rest of this study is arranged as follows. Section 2 gives detailed definition of the model. e expression of channel impulse response (CIR) is given in Section 3. In Section 4, the statistical properties are calculated. Section 5 is the simulation results and the comparison with measurement data. e conclusion is presented in Section 6.
Notation: the superscript (·) T denotes the transpose operator. e superscript (·) * denotes the conjugate operator. |x| is the norm of a vector x. E · { } denotes the expectation operator. p(α, β) is the probability density function (PDF) of the VMF distribution with parameters α and β.

Nonstationary 3D RS-GBSM for V2V
MIMO Channels e presented 3D RS-GBSM is for MIMO V2V communication systems, and the system structure is shown in Figures 1 and 2. We assume that both of Tx and Rx are surrounded by several vehicles. For convenience on building the coordinates, we assume that Tx follows Rx, and it can be very easily adapted to other cases. e coordinates are as shown in Figure 2. According to the different positions of Tx and Rx on the slope, it can be divided into 2 scenarios: (1) all vehicles are on the slope; (2) Tx and surrounding vehicles are off the slope and moving towards the bottom end of the slope while Rx and surrounding vehicles are on the slope.

Model Specific Description.
e distribution of scatterers in the multicluster form at the Tx side is shown in Figure 3. Since the distribution of scatterers at the Rx side is the same as the Tx side, we only illustrate the Tx side. Taking cluster S 1 as an example, the scatterers are distributed around the main path L 1 and obeying the VMF distribution. Since vehicles have a certain volume, the radii of different scatterers to Tx and Rx are different. Limited by the radii, they are distributed between s (1) 1 and s (2) 1 , i.e., all scatterers are distributed within the outermost and innermost two spheres. For clearer description of the structure, the concerned paths and parameters are shown in Figure 1. Major parameters in this study are listed in Table 1. According to the realistic situation, both the LoS and non-LoS (NLoS) components are included. A two-sphere model is developed for the NLoS components. As shown in Figures 1 and 3, effective scatterers are distributed between two spheres, and they are distributed in the form of clusters. Taking the fourcluster scenario as an example, there are four clusters of scatterers at both ends of Tx and Rx, representing, respectively, four vehicles around Tx and Rx, and scatterers are distributed on the surface and in these vehicles. As mentioned above, each effective scatterer has a different radius, and the specific way of modeling is shown in Figure 3. In Figure 1, we assume that both of Tx and Rx are equipped with multiple antennas. e antenna elements are omnidirectional, and the numbers of antennas at Tx and Rx are denoted as n t and n r , respectively. In this model, taking 2 × 2 MIMO as an example, the two antennas of the Tx are labeled p 1 and p 2 , and the two antennas of the Rx are labeled l 1 and l 2 . In addition, we consider that both Tx and Rx are moving at high speeds, and the directions of the movements of the Tx and Rx are jointly defined by elevation angles and azimuth  International Journal of Antennas and Propagation 3 angles. e azimuth and elevation angles of the movement directions are recorded as c T (t), c R (t), ϕ T (t), and ϕ R (t), and the velocities of the movements are recorded with v T (t)and v R (t). As mentioned before, the effective scatterers are distributed over a section of the sphere and are distributed in the form of clusters, denoted by S n ij i (i � 1, . . . , N c and j � 1, 2), where i represents that the scatterers are located in the i th cluster, j � 1 represents that the scatterers are on the Tx sphere, and j � 2 represents that the scatterers are on the Rx sphere. In the i th cluster, n ij represents that the scatterer is the n ij th scatterer of this cluster. α n ij T (t) and β n ij T (t) are the azimuth AoD (AAoD) and elevation AoD (EAoD) of the path. α n ij R (t) and β n ij R (t) are the azimuth AoA (AAoA) and elevation AoA (EAoA) of the path. e scattering paths can be divided into three types: (1) rays that are only reflected once at the Tx end, and these paths are called L SB 1 . (2) Rays that are only reflected once at the Rx end, and these paths are called L SB 2 . (3) Rays that are reflected twice at both of the Tx and Rx ends, and these paths are called L DB . With the distances between adjacent antennas δ T and δ R and the radii of the two spheres R T and R R , we have δ T , δ R ≪ R T , R R .

Vehicle Trajectory.
We have considered not only the situation where vehicles move in a straight line but also the situation where one vehicle changes lane or overtakes. Vehicle's trajectory is shown in Figure 4. e initial time and final time are recorded as t ini and t fin . e location, velocity, and acceleration of the vehicle in the X-axis and Y-axis directions at the initial time and the final time are recorded as ( and € x denote distance, velocity, and acceleration, respectively.
We define the location of the vehicle in the X-axis and the Y-axis direction as and we define the time parameter matrix: fin t fin 1 5t 4 fin 4t 3 fin 3t 2 fin t fin 1 0 20t 3 fin 12t 2 fin 6t fin 2 0 0 Now, we obtain e elevation angles of movement direction of Tx and Rx rad e azimuth angles of movement direction of Tx and Rx rad e maximum Doppler shifts caused by Tx and Rx Hz e Doppler shifts caused by Tx and Rx Hz e horizontal orientation of antenna array of Tx and Rx rad α  International Journal of Antennas and Propagation Finally, the solution of the inhomogeneous linear equations (3a) and (3b) can be derived through the boundary conditions. Using this method, we can quickly obtain the track of vehicle's overtaking or lane-changing. However, the trajectory does not consider the dynamic characteristics of the vehicle, and the generated acceleration should be limited by the condition of vehicle engine, tire, ground friction, and other restrictions. Commonly, we have − 2.5 m/s 2 < € x < 2.5 m/s 2 and 2 m/s 2 < € y < 2 m/s 2 .

Channel Impulse Response
In this section, the detailed deduction of the CIR is given. Take path p 1 ⟶ l 1 as an example. Since scenario 1 and scenario 2 are very similar in derivation, we only give the derivation process for scenario 1. e CIR can be divided into three parts: LoS, single-bounced (SB), and double-bounced (DB) components. en, the CIR can be expressed as e Doppler shift caused by the Tx's motion is Correspondingly, the Doppler shift caused by the Rx's motion can be expressed as are the maximum Doppler shifts at time t, λ is the carrier wavelength, and v T (t) and v R (t) are the velocity vectors of the Tx and Rx. en, the three parts of the CIR can be rewritten as where N i1 is the number of effective scatterers in the i th cluster on the Tx sphere, N i 2 is the number of effective scatterers in the i th cluster on the Rx sphere, K is the Ricean factor, and SB 1 , SB 2 , and DB, respectively, represent L SB 1 , L SB 2 , and L DB which have been introduced in Section 2. Parameters ε SB i and ε DB are related to the power of different paths, they control the proportions of power among paths L SB 1 , L SB 2 , and L DB , and they satisfy ε 2 SB 1 + ε 2 SB 2 + ε 2 DB � 1. ψ n ki and ψ n i 1 ,n k 2 are the random phases, and they are uniformly distributed within [− π, π). f LoS/SB i /DB T/R Doppler (t), respectively, represent the Doppler shifts at Tx and Rx of the LoS, SB, and DB components. e expressions are derived as (5a) and (5b). For In terms of SB and DB components, we have α T/R (t) � α We consider that the locations of Tx and Rx are known through the global positioning system; the location of p 1 is [x p 1 (t), y p 1 (t), z p 1 (t)], and the location of l 1 is [x l 1 (t), y l 1 (t), z l 1 (t)]. In scenario 1, Tx, Rx, and surrounding vehicles are constantly moving on the slope. At initial time, x p 1 (t) � 0, y p 1 (t) � δ T cos θ T (t)/2, z p 1 (t) � δ T sin θ T (t)/2, x l 1 � D(t)cos θ slope , y l 1 (t) � δ R cos θ R (t)/2, and z l 1 (t) � δ R sin θ R (t)sin θ slope /2. e location of S (n i 1 ) i can be denoted as [x s 1 (t), y s 1 (t), z s 1 (t)]: International Journal of Antennas and Propagation 5 en, we can easily calculate d(p 1 , l 1 )(t), d(p 1 , S In DB paths, EAoD/AAoD and EAoA/AAoA are independent, whereas in SB 1 and SB 2 paths, they are related. e generation of angles is introduced in Section 4.

Space-Time Correlation Function.
e normalized space-time correlation function (STCF) between path p 1 ⟶ l 1 and p 2 ⟶ l 2 is defined as e STCF can be written into (a) STCF of LoS component: we also only give the derivation for scenario 1 because it is similar to the scenario 2. e STCF of the LoS component is where d(p 1 , l 1 )(t) and d(p 2 , l 2 )(t) can be calculated as (9). f LoS T Doppler (t) and f LoS R Doppler (t) have been expressed in Section 3. (b) STCF of SB components: we consider that the EAoD and the AAoD in SB 1 and EAoA and AAoA in SB 2 are contributed to VMF contribution; thus, the general PDF of these angles is p(α, β) � k cos β 0 4π sinh k × e k cos β 0 cos β cos α− α 0 ( )+sin β 0 sin β , (13) where α, β ∈ [ − π, π ), and α 0 , β 0 ∈ [ − π, π ). Since there are N c clusters of scatterers, there are also N c groups of α 0 and β 0 corresponding to p i (α, β) (i �1, . . ., N c ). Next, if the number of the scatterers tends to infinity, T/R can be replaced by continuous random variables α T/R and β T/R . en, the STCF of SB i component becomes where , and the calculations are similar to (9).
(c) STCF of DB components: where

Doppler Power Spectrum Density.
If we apply the Fourier transform to the space-time correlation function, then we can get the corresponding Doppler power spectrum density (PSD).

Simulation Results and Analysis
In this section, several statistical characteristics of the proposed channel model are simulated and analyzed, including the spatial ACF, temporal ACF, and Doppler PSD. Based on measurement data in [28,52] and measurement methods in [53][54][55], the simulation parameters at the initial time are set and listed in Table 2. As for the scenario with slope scenario 2 is chosen in simulations. Due to the energy of the scattering path is much smaller than that contained in the LoS path, the LoS path is ignored in the spatial and temporal ACFs simulation. Besides, since the above measurement campaigns are completed on common roads without slope, in CIR, timevarying ACF, and Doppler PSD simulations, the propagation scenario is reduced to a level load.

Simulation of CIR.
As mentioned above, in the realistic scenario, the velocity and position vectors of the vehicle are constantly changing; thus, the nonstationary model is closer to reality. e tapped delay line (TDL) model is adopted in the CIR simulation. In Figure 5, it can be clearly seen that the locations of different clusters change with time. It needs to be mentioned that the reality scenario is more complicated, since the power of the ray is constantly changing, and new scatterers will appear that accompany existed scatterers' disappearance. In Figure 6, we give a snapshot of the TDL at t � 5 s. It can be seen that due to the setting of Tx, Rx, and other vehicles' velocities and locations, different clusters gather together after a certain period of time.

Spatial ACF.
is part is the simulation result of spatial ACF. Figure 7 shows the spatial ACFs of scenarios with and without slope at Tx end. e spacing between the two receiving antennas at Rx is fixed, and the antenna spacing in Figure 7 refers to the spacing between adjacent transmitting antennas at Tx. For example, in our simulation, there are N c � 4 clusters of scatterers around the Tx and Rx, respectively, and each clusters contains 10 effective scatterers (N i 1 � N j 2 � 10). Since there are at total N c × N i 1 × N c × N j 2 DB paths, the spatial diversity is huge. As a result, the spatial ACF of the channel is therefore relatively low. Figure 7 presents normalized spatial ACFs for scenarios including slope with different tilt angles, i.e., θ slope � 0, θ slope � π/36, θ slope � π/6, and θ slope � π/3. e figure shows that the tilt angle of the slope has little effect on the spatial ACF. As the tilt angle of slope increases, the main lobe of the spatial ACF becomes narrower.

e ACF at Initial Time.
When simulating the temporal ACF, it is necessary to, respectively, fix the spacing of adjacent transmitting and receiving antennas at Tx and Rx. We take δ T � δ R � λ in the simulations. e generation of International Journal of Antennas and Propagation the scatterers is the same as that in the simulation of the spatial ACF. e specific results are shown in Figure 8.
In the VMF distribution, the number of α 0 is 2 × N c because there are N c clusters of scatters near Tx and another N c clusters of scatters near Rx, the values of them are listed in Table 2, and the values of β 0 are set to zero. Figure 8 shows the temporal ACFs for scenarios with a slope of different tilt angles, and the angles are the same as above. It can be seen from the figure that as the tilt angle of slope increases, the main lobe of the temporal ACF becomes wider.

Time-Varying ACF.
In this simulation, we assume that the speeds of the Tx, Rx, and vehicles around them are 30 m/s, 40 m/s, and 35 m/s, respectively. e results of the simulation are shown in Figures 9 and 10.
From Figure 9, we can observe the nonstationary feature of the channel. It also verifies that when the speed directions of Tx and Rx coincide, the channel can be regarded as a stationary channel in a short interval (within 1 second). Figure 10 illustrates more detailed time-variant features of the ACFs at different moments, e.g., at the moments t � 0 s, 1 s, 3 s, and 5 s. In the simulated scenario set, different vehicles are moving away from each other as time increases. erefore, the different scattering paths are constantly moving closer together, as shown in Figure 11, the AAoA of SB 1 becomes closer to the LoS path (0 radians). And the difference between the angles of the scattering paths and the delays are continuously decreasing. Eventually, it leads to an increase in the main lobe of the ACF.

Doppler PSD.
e simulated Doppler PSDs are shown in Figures 12 and 13. e result of Figure 12 is consistent with Figure 5. As time increases, Tx and Rx keep moving away from each other, and different clusters gather with each other. It results in a decrease in the width of Doppler frequency.
We have N i 1 � N j 2 � 10, so there should be a total of N clusters of impulses in Figure 13, i.e., N � N 1 + N 2 + N 3 , Figure 13, only N 1 + N 2 clusters of SB path impulses can be seen because the energy of DB path is too low. Furthermore, in each cluster of scatterers, all scatterers are surrounding the same direction, so as it can be seen in Figure 13 that these impulses are in the form of clusters. Besides, the Doppler PSD provides help to verify the correctness of the model. After comparison and verification, the Doppler PSDs shown in Figures 12 and 13 are consistent with the theoretical analysis.

Comparison with Measurements.
In this section, we compare the simulation results of the model with measurement data to verify the accuracy of the model. We choose the parameter and fitting cure given in [52] and adjust the vehicle's distribution and movement in the model to match the measurement scenario. After adjusting the model, we simulate RMS delay spread, RMS Doppler spread, and RMSE to be compared with the measurement data.

Adjustments of Scenario.
e adjusted scenario is shown in Figure 14. Since the measurement environments are level roads, we set the tilt angle of slope to zero. At the initial time, Tx and Rx are moving on two adjacent lanes, then Rx accelerates and overtakes Tx. e velocity of different vehicles varies according to different scenarios.
In [52], measurement campaigns are in four different scenarios, highway, urban, suburb, and municipal lake. In the lake scenario, there is very little scattering caused by the environment because there are almost no buildings in this scenario, and the small-scale fading of the channel mainly comes from the scattering caused by the surrounding vehicles. In other measurement scenarios, the environment causes severe channel fading. erefore, we added  D 2 ) (100, 100) m (R T , R R ) ( (15,20), (15,20)) m ( (15,20), (15,20)    Besides, it can be divided into overtaking (OT) and nonovertaking (NOT) cases according to the distance between Tx and Rx. In the case of OT, the distance between Tx and Rx is very close, less than 11 m. While in the NOT case, they keep a certain distance.

RMS Delay Spread.
In wireless communication, RMS delay spread is a statistical parameter that can describe the delay characteristics of the radio channel. It is often used in the evaluation of channel characteristics. In V2V communication, the multipath effect caused by the scattering between vehicles and scattering from the environment caused small-scale fading. RMS delay spread can reflect this multipath effect. e complex time-varying CIR can be derived by an inverse Fourier transform of the channel transfer function (6a)-(6c). en, we can calculate the time-variant power delay profile (PDP) by averaging the magnitude squared of the CIR.
where a i is the complex coefficient of each delay path, and τ i is the excess delay of the i th path. en, we can get the RMS delay T rms as follows: where P A (τ) is the averaged power delay profile. e comparison results are shown in Figure 15. e fitting curves corresponding to different scenarios are given   in [52]. We put simulation results, measurement data, and fitting curves together for comparison. In the following section, RMSE will be used for evaluation.

RMS Doppler Spread.
In V2V communication, vehicles are moving at high speed, and the speed direction of vehicles may change, so it causes a strong Doppler effect. erefore, RMS Doppler spread is also applied to the description of V2V channel delay characteristics.
We can get the RMS Doppler spread D rms as where P Bm represents the integration of scattering function.
e comparison results are shown in Figure 16. e comparison of simulation results, measurement data, and fitting curves are given in the figure.

RMSE.
To verify that the model can fit the measurement data better than the fitting curve, we calculate the RMSE of delay and Doppler spread.
where y i denotes the measurement data, and x i denotes the simulation result or the fitting curve. RMSE can reflect the fitting degree of simulation results and measurement data. A smaller RMSE means a better fitting. e calculation results are shown in Table 3. When the difference of RMSE is negative, it indicates that the model can better fit the measurement data than the fitting curve. As can be seen in the table, the model has a smaller RMSE in most cases, and the model has a similar degree to the fitting curve in the remaining scenarios.

Conclusion
In this study, we developed an analytical multicluster nonstationary 3D MIMO channel model for the V2V communication system, which is suitable for the scenario where there are a certain number of vehicles near Tx and Rx on or near a slope. After simplification, the model can be applied to the normal road without slope. In addition, we also considered the impact of vehicle changing lanes. Furthermore, we derived and simulated statistical properties of the proposed channel such as the spatial ACF, temporal ACF, time-varying ACF, and Doppler PSD. Finally, we compared the simulation results with the measurement data in four scenarios. e comparison results show that the proposed model can mimic the measurement data better than the fitting curve. ese simulation results have shown the utility of the proposed model.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that they have no conflicts of interest.