# A 3D Non-Stationary Channel Model with Moving Mobile Station in Rectangular Tunnel

Authors: Zhicheng Ren; Fangqi Zhang; Guoxin Zheng; Asad Saleem; Ke Guan
Journal: International Journal of Antennas and Propagation
 (2019)

---

## Abstract

This paper presents a three-dimensional (3D) geometry-based stochastic model (GBSM) for capturing the non-stationarity of radio channel at 1.8 GHz in a rectangular tunnel. A time-variant (TV) complex channel gain is derived for obtaining the statistical properties in time, frequency, and spatial domains such as the time-variant autocorrelation function (TV-ACF), the time-variant Doppler power spectral density (TV-DPSD), and the time-variant spatial cross-correlation function (TV-CCF), respectively. Then the TV channel statistical properties at different time instants can be extracted and the non-stationary channel characteristics caused by the TV scattering environment are thoroughly discussed. Furthermore, three cases including “approach”, “arrival”, and “away” are set to allow a comprehensive study on how the DPSD behaves with the relative position between transmitter and receiver. The reliability of proposed 3D GBSM is highlighted by a good agreement with the measured result in terms of the correlation function.


---

## Body

## 1. Introduction

For realizing a more stable and reliable train-ground communication, the multiple-input multiple-output (MIMO) technology working at 1.8 GHz has been gradually utilized in metro communication system. The MIMO channel model is thus a prerequisite for predicting the coverage of wireless link and optimizing the system performance. The geometry-based stochastic model (GBSM) [1–3] which combines the geometry of the environment and the statistical information is regarded as a flexible way to describe the MIMO channel characteristics. Considering the propagation in the elevation plane, lots of attention has been paid to three-dimensional (3D) channel models [4–6]. In [4], a concentric-cylinder model was derived for wideband MIMO mobile-to-mobile channel, and the independent azimuth angle and elevation angle was analyzed separately. To study the joint effect of the azimuth and elevation angle, 3D two-sphere models have been proposed in [5, 6]. Several GBSM based standard channel models also have been proposed, e.g., 3GPP-3D model, IMT-Advanced model, and WINNER+ model. These standard models were built with the angle parameter obeying Laplacian, Gaussian, and Wrapped Gaussian distributions, respectively. Furthermore, these models were evolved to support to non-stationary scenario by replacing the model parameters with the time-variant (TV) ones [7, 8]. In spite of this, these standard models are still focused on propagation scenarios like Indoor, Urban, Rural, and Suburban; mobile propagation scenarios are not involved in these models. Another modeling approach in [9, 10] employed random mobility trajectories of the mobile station (MS) to inject the non-stationarity into the characteristics of the channel, but without the consideration of the regular movement in confined space, especially in tunnel scenario. It is worth mentioning that the radio signal propagations inside tunnel will suffer more reflections and scattering due to the confined space, rough interior tunnel walls, and waveguide effect [11]. To acquire better performance of metro communication system, it is essential to build accurate tunnel channel models to comprehensively understand the statistical properties of tunnel channels.Few studies were devoted to the non-stationary channel models for the tunnel environment. The regular-shaped GBSM in cutting and irregular-shaped GBSM in suburb and urban scenarios were proposed in [12]. In [13], a geometry-based single-bounce model of MIMO channel in the tunnel has been proposed in order to investigate the factors which influence the channel correlation, but it is estimated without considering the elevation angles. A geometry-based MIMO car-to-car channel model is applied for a semi-circular tunnel [14]. However, this model only considered the local scatterers between two cars which both move at the same speed. The local angle of arrival (AOA) and angle of departure (AOD) were thus time-invariant, and the non-stationarity was unable to be characterized. The two-dimensional (2D) high-speed train model in [15] considered the non-stationary characteristic caused by the variation of the velocity and the angle of motion (AOM), but the temporal variations of AOA which contribute to the non-stationary characteristic were neglected.However, the modeling of AOA and AOD plays an import role in the non-stationary mobility model. Empirical stochastic distributions, e.g., the uniform distribution [16, 17], Von Mises distribution [18, 19], and Von Mises-Fisher distribution [20], have always been assumed to characterize the angular parameter to achieve the closed-form representation of channel statistical properties. But the AOAs in [16–20] are not modeled as time-variant functions. It is not suitable for the real-time mobile behavior modeling; the temporal variations of the angle parameters need to be included in these models.Inspired by the issues mentioned before, we present a 3D geometrical model to describe the non-stationary channel caused by the movements of the MS in the tunnel environment. By introducing the TV AOA and propagation length, we derive the TV channel gain based on the single-bounce scattering. Furthermore, in the non-stationary process, the phase of the complex channel gain is no longer rotated with the TV Doppler frequency multiplied by time, but the integral of the TV Doppler frequency. The former multiplicative relationship can only be obtained by developing the TV phase in a second-order Taylor series aroundt = 0 [21]. By employing the analytical method of the non-stationary processes, the expressions and analytical results are given for the time-variant autocorrelation function (TV-ACF), the time-variant Doppler power spectral density (TV-DPSD), and the time-variant spatial cross-correlation function (TV-CCF). Moreover, the correctness of the derivations is verified by the simulation results. Finally, a measurement is conducted in a subway-like tunnel and the measured results are used to validate the effectiveness of the proposed model.The paper is organized as follows. Section2 proposes a new 3D non-stationary GBSM, and the statistical properties of the model are derived. The theoretical results are depicted and compared with the corresponding simulation in Section 3. In Section 4, the measurement results are discussed with the proposed model. Finally, Section 5 concludes the paper and outlines future research lines.


## 2. The 3D Non-Stationary GBSM

A typical propagation scenario for the metro environment is shown in Figure1. The downlink is from the base station (BS) to the moving MS. The symbols of BST and MSR denote the BS transmitter (Tx) and the MS receiver (Rx), respectively. The proposed model generates the empty straight tunnel scattering environment with a length of L. The tunnel cross section is rectangular with a width of W and a height of H. Figure 2 shows the distribution of scatterers on the inner surfaces of the tunnel wall. In order to study the general cases of the scattering inside the tunnel with the moving mobile station, it is assumed that the material of the tunnel wall is uniform and the scatterer distribution follows uniform distribution.Figure 1
Propagation environment inside tunnels.Figure 2
The distribution of scatterers.The 3D MIMO non-stationary GBSM is described as follows. The measurement inside the tunnel [22] demonstrates that most of the propagation wave energy was from the line-of-sight (LOS) and single-bounced (SB) components. Therefore, the LOS and SB components are under consideration in our proposed model. The signals received at the moving Rx are reflected from the infinite scatterers which are randomly distributed along the walls, ceiling, and ground of tunnel. For the simplified description, (xm, yn, zo) are used to describe the positions of the scatterers S(m,n,o) on the inner surfaces of the tunnel wall for m=1,2,…,M, n=1,2,…,N, and o=1,2,…,O. Since the upper limitM,N,O with each tending to infinity, three discrete random variables xm, yn, and zo can be denoted as the continuous variablesx,y, and z to acquire the theoretical model.


### 2.1. Channel Gain

The channel gain between thep-th Tx element BSpT and theq-th Rx element MSqR is demonstrated as the superposition of the LOS and the SB components, and it can be expressed as(1)hpqt=hpqLOSt+hpqSBtFor the SB component, the channel gainhpqSBt can be formulated as the summation of several multipath components from M×N×O scatterers as given below(2)hpqSBt=ΩpqKpq+1limM→∞N→∞O→∞⁡1MNO∑m,n,o=1M,N,Oe-j2π/λDpqm,n,otej2π∫0tfRSBudu+θm,n,owhere(3)Dpqm,n,ot=DpTm,n,ot+DqRm,n,ot(4)DpTm,n,ot=xT-xm2+yT-yn2+zT-zo2-δT2MT-2p+1cosϕT·cos⁡βTm,n,ot·cos⁡γT-αTm,nt+sinϕT·sin⁡βTm,n,ot(5)DqRm,n,ot=xR+vt-xm2+yR-yn2+zR-zo2-δR2MR-2q+1cosϕR·cos⁡βRm,n,ot·cos⁡γR-αRm,nt+sinϕR·sin⁡βRm,n,ot(6)fRSBt=fRmaxcos⁡αRm,nt-φvRcos⁡βRm,n,otThe key parameters in equations are defined in Table 1. The symbol Kpq in (2) represents the Rice factor, which is defined as Kpq=hpqLOSt2/EhpqSBt2. Ωpq denotes the mean power between thep-th Tx element and theq-th Rx element. In (4), (xT,yT,zT) and (xR,yR,zR) are used to represent the positions of the BST and MSR. Using the trigonometric relationship, the azimuth angle of arrival and elevation angle of arrival can be expressed as(7)
α
R
m
,
n
t
=
arccos
⁡
x
m
-
x
R
-
v
t
x
m
-
x
R
-
v
t
2
+
y
n
-
y
R
2
,
y
n
≥
y
R
-
arccos
⁡
x
m
-
x
R
-
v
t
x
m
-
x
R
-
v
t
2
+
y
n
-
y
R
2
,
y
n
<
y
R
(8)
β
R
m
,
n
,
o
t
=
arccos
⁡
x
m
-
x
R
-
v
t
2
+
y
n
-
y
R
2
x
m
-
x
R
-
v
t
2
+
y
n
-
y
R
2
+
z
o
-
z
R
2
,
z
o
≥
z
R
-
arccos
⁡
x
m
-
x
R
-
v
t
2
+
y
n
-
y
R
2
x
m
-
x
R
-
v
t
2
+
y
n
-
y
R
2
+
z
o
-
z
R
2
,
z
o
<
z
RThe azimuth angle of departure αTm,nt and the elevation angle of departure βTm,n,ot can be derived in the similar way. In (6), fRmax=v/λ indicates the maximum Doppler shift and fRSBt denotes the TV Doppler shift of the SB component caused by the movement of the Rx, where λ is the wavelength.Table 1
Definition of parameters.
Parameters Definition D p T m , n , o t , D q R m , n , o t the SB path distanced(BSpT,S(m,n,o)) and d(Sm,n,o,MSqR) D p q t the LOS path distanced(BSpT,MSqR) α T m , n, βTm,n,ot the azimuth and elevation angle of departure α R m , n, βRm,n,ot the azimuth and elevation angle of arrival δ T,δR the antenna element spacings of the Tx and Rx M T, MR the antenna element number of the Tx and Rx γ T, γR the azimuth angle of the Tx and Rx ϕ T, ϕR the elevation angle of the Tx and Rx v, φvR the moving speed and the moving angle of the RxThe LOS componenthpqLOSt of the total channel gain in (1) can be written as (9)hpqLOSt=ΩpqKpqKpq+1e-j2π/λDpqtej2π∫0tfRLOSudu+θLOSwhere(10)Dpqt=xT-xR-vt2+yT-yR2+zT-zR2-δT2MT-2p+1cosϕTcosγT-δR2MR-2q+1cosϕRcos⁡αTLOSt-γR(11)fRLOSt=fRmaxcos⁡αRLOSt-φvRcosβRLOStIt is worth mentioning that the θm,n,o in (2) and the θLOS in (9) are the initial phases of the SB and the LOS component; each is assumed with uniform distribution over the interval (0,2π]. The symbol fRLOSt in (11) denotes the TV Doppler shift of the LOS component caused by the movement of the Rx. The TV angle parameters αTLOSt, αRLOSt, and βRLOSt also can be computed in the similar way as (7) and (8).


### 2.2. Statistical Properties of the Model

The typical TV channel statistical properties can be derived from the channel gain, such as the TV-ACF, TV-DPSD, and TV-CCF.


#### 2.2.1. TV-ACF

The TV-ACF of the subchannel from thep-th Tx element to theq-th Rx element is defined by(12)rpqt,τ=Ehpqthpq∗t-τΩpq=rpqLOSt,τ+rpqSBt,τwhere ∗ denotes the complex conjugate operator and E{·} is the statistical expectation operator.Substituting (2) and (9) into (12) and assuming that the propagation distance is invariant for short time difference τ and the phase rotated with the integral of the TV Doppler frequency has a prominent impact on the TV-ACF, the TV-ACF of the LOS and the SB component can be obtained as follows:(13)rpqLOSt,τ=KpqKpq+1ej2π∫t-τtfRmaxcos⁡αRLOSu-φvRcosβRLOSudu(14)rpqSBt,τ=1Kpq+1∫0L∫-W/2W/2∫0Hfxyzx,y,zej2π∫t-τtfRmaxcos⁡αRx,yu-φvRcosβRx,y,zududxdydzwhere fxyzx,y,z is the joint probability density function (PDF) of the continuous random valuesx,y, andz.


#### 2.2.2. TV-DPSD

For non-stationary channels, the TV-DPSD can be obtained from the TV-ACF by applying the Fourier transform with regard to time differenceτ, as given below(15)Spqt,f=∫-∞∞rpqt,τe-j2πfτdτIt can be executed to reflect the variations in the power spectrum along the Doppler frequency throughout the movement of the Rx.


#### 2.2.3. TV-CCF

The normalized TV-CCF between two arbitrary channels in the spatial domain is(16)ρpq,p′q′t,δT,δR=Ehpqthp′q′∗tΩpqΩp′q′=ρpq,p′q′LOSt,δT,δR+ρpq,p′q′SBt,δT,δRFor the case of the LOS component, it can be defined as(17)ρpq,p′q′LOSt,δT,δR=KpqKpq+1e-j2π/λδTp-p′cosϕTcosγT·e-j2π/λδRq-q′cosϕRcos⁡αTLOSt-γRFor the case of the SB component, it can be defined as(18)ρpq,p′q′SBt,δT,δR=1Kpq+1∫0L∫-W/2W/2∫0Hfxyzx,y,z·e-j2π/λδTp-p′cosϕTcosβTx,y,ztcosγT-αTx,yt·e-j2π/λδTp-p′sinϕTsinβTx,y,zt·e-j2π/λδRq-q′cosϕRcosβRx,y,ztcosγR-αRx,yt·e-j2π/λδRq-q′sinϕRsinβRx,y,ztdxdydzMoreover, the CCF at BST and MSR can be denoted as ρpq,p′q′SBt,δT and ρpq,p′q′SBt,δR, respectively. As an example, by substituting δT=0 in (18), ρpq,p′q′SBt,δR can be obtained as follows:(19)ρpq,p′q′SBt,δR=1Kpq+1∫0L∫-W/2W/2∫0Hfxyzx,y,z·e-j2π/λδRq-q′cosϕRcosβRx,y,ztcosγR-αRx,yt·e-j2π/λδRq-q′sinϕRsinβRx,y,ztdxdydzThe uniform distribution is assumed in scattering model; therefore the joint PDF can be computed as [23](20)fxyzx,y,z=2LW2+4z2arctan2H/W


## 2.1. Channel Gain

The channel gain between thep-th Tx element BSpT and theq-th Rx element MSqR is demonstrated as the superposition of the LOS and the SB components, and it can be expressed as(1)hpqt=hpqLOSt+hpqSBtFor the SB component, the channel gainhpqSBt can be formulated as the summation of several multipath components from M×N×O scatterers as given below(2)hpqSBt=ΩpqKpq+1limM→∞N→∞O→∞⁡1MNO∑m,n,o=1M,N,Oe-j2π/λDpqm,n,otej2π∫0tfRSBudu+θm,n,owhere(3)Dpqm,n,ot=DpTm,n,ot+DqRm,n,ot(4)DpTm,n,ot=xT-xm2+yT-yn2+zT-zo2-δT2MT-2p+1cosϕT·cos⁡βTm,n,ot·cos⁡γT-αTm,nt+sinϕT·sin⁡βTm,n,ot(5)DqRm,n,ot=xR+vt-xm2+yR-yn2+zR-zo2-δR2MR-2q+1cosϕR·cos⁡βRm,n,ot·cos⁡γR-αRm,nt+sinϕR·sin⁡βRm,n,ot(6)fRSBt=fRmaxcos⁡αRm,nt-φvRcos⁡βRm,n,otThe key parameters in equations are defined in Table 1. The symbol Kpq in (2) represents the Rice factor, which is defined as Kpq=hpqLOSt2/EhpqSBt2. Ωpq denotes the mean power between thep-th Tx element and theq-th Rx element. In (4), (xT,yT,zT) and (xR,yR,zR) are used to represent the positions of the BST and MSR. Using the trigonometric relationship, the azimuth angle of arrival and elevation angle of arrival can be expressed as(7)
α
R
m
,
n
t
=
arccos
⁡
x
m
-
x
R
-
v
t
x
m
-
x
R
-
v
t
2
+
y
n
-
y
R
2
,
y
n
≥
y
R
-
arccos
⁡
x
m
-
x
R
-
v
t
x
m
-
x
R
-
v
t
2
+
y
n
-
y
R
2
,
y
n
<
y
R
(8)
β
R
m
,
n
,
o
t
=
arccos
⁡
x
m
-
x
R
-
v
t
2
+
y
n
-
y
R
2
x
m
-
x
R
-
v
t
2
+
y
n
-
y
R
2
+
z
o
-
z
R
2
,
z
o
≥
z
R
-
arccos
⁡
x
m
-
x
R
-
v
t
2
+
y
n
-
y
R
2
x
m
-
x
R
-
v
t
2
+
y
n
-
y
R
2
+
z
o
-
z
R
2
,
z
o
<
z
RThe azimuth angle of departure αTm,nt and the elevation angle of departure βTm,n,ot can be derived in the similar way. In (6), fRmax=v/λ indicates the maximum Doppler shift and fRSBt denotes the TV Doppler shift of the SB component caused by the movement of the Rx, where λ is the wavelength.Table 1
Definition of parameters.
Parameters Definition D p T m , n , o t , D q R m , n , o t the SB path distanced(BSpT,S(m,n,o)) and d(Sm,n,o,MSqR) D p q t the LOS path distanced(BSpT,MSqR) α T m , n, βTm,n,ot the azimuth and elevation angle of departure α R m , n, βRm,n,ot the azimuth and elevation angle of arrival δ T,δR the antenna element spacings of the Tx and Rx M T, MR the antenna element number of the Tx and Rx γ T, γR the azimuth angle of the Tx and Rx ϕ T, ϕR the elevation angle of the Tx and Rx v, φvR the moving speed and the moving angle of the RxThe LOS componenthpqLOSt of the total channel gain in (1) can be written as (9)hpqLOSt=ΩpqKpqKpq+1e-j2π/λDpqtej2π∫0tfRLOSudu+θLOSwhere(10)Dpqt=xT-xR-vt2+yT-yR2+zT-zR2-δT2MT-2p+1cosϕTcosγT-δR2MR-2q+1cosϕRcos⁡αTLOSt-γR(11)fRLOSt=fRmaxcos⁡αRLOSt-φvRcosβRLOStIt is worth mentioning that the θm,n,o in (2) and the θLOS in (9) are the initial phases of the SB and the LOS component; each is assumed with uniform distribution over the interval (0,2π]. The symbol fRLOSt in (11) denotes the TV Doppler shift of the LOS component caused by the movement of the Rx. The TV angle parameters αTLOSt, αRLOSt, and βRLOSt also can be computed in the similar way as (7) and (8).


## 2.2. Statistical Properties of the Model

The typical TV channel statistical properties can be derived from the channel gain, such as the TV-ACF, TV-DPSD, and TV-CCF.


### 2.2.1. TV-ACF

The TV-ACF of the subchannel from thep-th Tx element to theq-th Rx element is defined by(12)rpqt,τ=Ehpqthpq∗t-τΩpq=rpqLOSt,τ+rpqSBt,τwhere ∗ denotes the complex conjugate operator and E{·} is the statistical expectation operator.Substituting (2) and (9) into (12) and assuming that the propagation distance is invariant for short time difference τ and the phase rotated with the integral of the TV Doppler frequency has a prominent impact on the TV-ACF, the TV-ACF of the LOS and the SB component can be obtained as follows:(13)rpqLOSt,τ=KpqKpq+1ej2π∫t-τtfRmaxcos⁡αRLOSu-φvRcosβRLOSudu(14)rpqSBt,τ=1Kpq+1∫0L∫-W/2W/2∫0Hfxyzx,y,zej2π∫t-τtfRmaxcos⁡αRx,yu-φvRcosβRx,y,zududxdydzwhere fxyzx,y,z is the joint probability density function (PDF) of the continuous random valuesx,y, andz.


### 2.2.2. TV-DPSD

For non-stationary channels, the TV-DPSD can be obtained from the TV-ACF by applying the Fourier transform with regard to time differenceτ, as given below(15)Spqt,f=∫-∞∞rpqt,τe-j2πfτdτIt can be executed to reflect the variations in the power spectrum along the Doppler frequency throughout the movement of the Rx.


### 2.2.3. TV-CCF

The normalized TV-CCF between two arbitrary channels in the spatial domain is(16)ρpq,p′q′t,δT,δR=Ehpqthp′q′∗tΩpqΩp′q′=ρpq,p′q′LOSt,δT,δR+ρpq,p′q′SBt,δT,δRFor the case of the LOS component, it can be defined as(17)ρpq,p′q′LOSt,δT,δR=KpqKpq+1e-j2π/λδTp-p′cosϕTcosγT·e-j2π/λδRq-q′cosϕRcos⁡αTLOSt-γRFor the case of the SB component, it can be defined as(18)ρpq,p′q′SBt,δT,δR=1Kpq+1∫0L∫-W/2W/2∫0Hfxyzx,y,z·e-j2π/λδTp-p′cosϕTcosβTx,y,ztcosγT-αTx,yt·e-j2π/λδTp-p′sinϕTsinβTx,y,zt·e-j2π/λδRq-q′cosϕRcosβRx,y,ztcosγR-αRx,yt·e-j2π/λδRq-q′sinϕRsinβRx,y,ztdxdydzMoreover, the CCF at BST and MSR can be denoted as ρpq,p′q′SBt,δT and ρpq,p′q′SBt,δR, respectively. As an example, by substituting δT=0 in (18), ρpq,p′q′SBt,δR can be obtained as follows:(19)ρpq,p′q′SBt,δR=1Kpq+1∫0L∫-W/2W/2∫0Hfxyzx,y,z·e-j2π/λδRq-q′cosϕRcosβRx,y,ztcosγR-αRx,yt·e-j2π/λδRq-q′sinϕRsinβRx,y,ztdxdydzThe uniform distribution is assumed in scattering model; therefore the joint PDF can be computed as [23](20)fxyzx,y,z=2LW2+4z2arctan2H/W


## 2.2.1. TV-ACF

The TV-ACF of the subchannel from thep-th Tx element to theq-th Rx element is defined by(12)rpqt,τ=Ehpqthpq∗t-τΩpq=rpqLOSt,τ+rpqSBt,τwhere ∗ denotes the complex conjugate operator and E{·} is the statistical expectation operator.Substituting (2) and (9) into (12) and assuming that the propagation distance is invariant for short time difference τ and the phase rotated with the integral of the TV Doppler frequency has a prominent impact on the TV-ACF, the TV-ACF of the LOS and the SB component can be obtained as follows:(13)rpqLOSt,τ=KpqKpq+1ej2π∫t-τtfRmaxcos⁡αRLOSu-φvRcosβRLOSudu(14)rpqSBt,τ=1Kpq+1∫0L∫-W/2W/2∫0Hfxyzx,y,zej2π∫t-τtfRmaxcos⁡αRx,yu-φvRcosβRx,y,zududxdydzwhere fxyzx,y,z is the joint probability density function (PDF) of the continuous random valuesx,y, andz.


## 2.2.2. TV-DPSD

For non-stationary channels, the TV-DPSD can be obtained from the TV-ACF by applying the Fourier transform with regard to time differenceτ, as given below(15)Spqt,f=∫-∞∞rpqt,τe-j2πfτdτIt can be executed to reflect the variations in the power spectrum along the Doppler frequency throughout the movement of the Rx.


## 2.2.3. TV-CCF

The normalized TV-CCF between two arbitrary channels in the spatial domain is(16)ρpq,p′q′t,δT,δR=Ehpqthp′q′∗tΩpqΩp′q′=ρpq,p′q′LOSt,δT,δR+ρpq,p′q′SBt,δT,δRFor the case of the LOS component, it can be defined as(17)ρpq,p′q′LOSt,δT,δR=KpqKpq+1e-j2π/λδTp-p′cosϕTcosγT·e-j2π/λδRq-q′cosϕRcos⁡αTLOSt-γRFor the case of the SB component, it can be defined as(18)ρpq,p′q′SBt,δT,δR=1Kpq+1∫0L∫-W/2W/2∫0Hfxyzx,y,z·e-j2π/λδTp-p′cosϕTcosβTx,y,ztcosγT-αTx,yt·e-j2π/λδTp-p′sinϕTsinβTx,y,zt·e-j2π/λδRq-q′cosϕRcosβRx,y,ztcosγR-αRx,yt·e-j2π/λδRq-q′sinϕRsinβRx,y,ztdxdydzMoreover, the CCF at BST and MSR can be denoted as ρpq,p′q′SBt,δT and ρpq,p′q′SBt,δR, respectively. As an example, by substituting δT=0 in (18), ρpq,p′q′SBt,δR can be obtained as follows:(19)ρpq,p′q′SBt,δR=1Kpq+1∫0L∫-W/2W/2∫0Hfxyzx,y,z·e-j2π/λδRq-q′cosϕRcosβRx,y,ztcosγR-αRx,yt·e-j2π/λδRq-q′sinϕRsinβRx,y,ztdxdydzThe uniform distribution is assumed in scattering model; therefore the joint PDF can be computed as [23](20)fxyzx,y,z=2LW2+4z2arctan2H/W


## 3. Simulation Results and Validation

Here, we investigate the statistical properties of the proposed model in detail. A 150 m long rectangular tunnel with the height of 5 m and the width of 5.2 m is considered.


### 3.1. Simulation Model

By assuming the infinite number of the scatterers (M,N,O→∞), the analytical model described in (2) is hard to realize. However, a corresponding simulation model with a finite number of scatterers can be developed. In terms of the model parameters, the discrete random values {xm,yn,zo} are generated by discrete uniform random parameters, while the other parameters are the same as the theoretical model, summarized in Table 2. To achieve a better approximation to the statistical properties of the theoretical model, the simulation results are obtained by averaging the results of multiple simulation runs (up to 104), and the number of scatterers per simulation run is selected as M×N×O=2000. Furthermore, the statistical properties of the simulation model are shown as (21)-(23). For TV-ACF, the simulation model is expressed as(21)r~pqSBt,τ=1Kpq+1·1MNO⁡∑m,n,o=1M,N,Oej2π∫t-τtfRmaxcos⁡αRm,nu-φvRcosβRm,n,ouduFor TV-DPSD, the simulation model is given as(22)S~pqSBt,f=∫-∞∞r~pqSBt,τe-j2πfτdτFor TV-CCF, i.e., the CCF at MSR, the simulation model is written as follows:(23)ρ~pq,p′q′SBt,δR=1Kpq+1·1MNO⁡∑m,n,o=1M,N,Oe-j2π/λδRq-q′cosϕRcosβRm,n,otcosγR-αRm,nt·e-j2π/λδRq-q′sinϕRsinβRm,n,otTable 2
Simulation parameters.
Parameters Value M T, MR 2, 2 f c (GHz) 1.8 v T, vR (km/h) 0, 80 (L,W,H) (m) (150, 5.2, 5) ( x T , y T , z T ) (m) (75, 2, 3) ( x R , y R , z R ) (m) (0, 0, 3) γ T, ϕT (rad) 0, 0 γ R, ϕR (rad) π / 2,  0 Ω p q (W) 1


### 3.2. TV-ACF

By using (12)-(14), the numerical results obtained for the TV-ACF rpqt,τ in NLOS scenario are shown in Figure 3. For clearly comparison purpose, the absolute values of the theoretical TV-ACF at different time instantst = 0 s, 1.5 s, 3 s (corresponding to xR = 0 m, 33.3 m, 66.6 m) are extracted as shown in Figure 4. Also, the corresponding simulation results can be obtained by (21). It can be seen that a good agreement between the theoretical and simulated results is achieved, which verifies the derivation of the proposed model. Furthermore, different variation trends of the ACF are resulted from the dynamic instantaneous phase which varies with the TV Doppler frequency. As the Rx gets closer to the Tx (fixed location, xT = 75 m), the channel autocorrelation behaves more sensitive to the time difference (see, for example, the curve att = 3 s). From Figure 5, it is obvious that the TV-ACF of the channel with the LOS component is higher than the case without the LOS component. In addition, the non-stationarity of channel characteristic is not notable in LOS dominated propagation condition.Figure 3
Absolute values of the theoretical TV-ACFrpqt,τ (Kpq = 0).Figure 4
Absolute values of the theoretical and simulated TV-ACFrpqt,τ at different time instants (Kpq = 0).Figure 5
Absolute values of the theoretical TV-ACFrpqt,τ with and without the LOS component at different time instants.


### 3.3. TV-DPSD

The theoretical TV-DPSD calculated by (15) is given in Figure 6. With the variation of the angle between the incoming wave and the direction of motion, the Doppler shift has a smooth drift from positive maximum to negative minimum. Figure 7 depicts three typical Rx locations with respect to Tx, the approach location (t = 0.45 s, xR = 10 m, Case I), the arrival location (t = 3.375 s, xR = 75 m, Case II), and the away location (t = 6.30 s, xR = 140 m, Case III). They are extracted from Figure 6 for further illustration as shown in Figure 8. A symmetric DPSD can be obtained in Case II while asymmetric DPSDs exist in Case I and Case III. In the non-isotropic environment like Case I and Case III, the AOA has a relatively concentrated distribution only in one direction. As a result, the energy focuses on the positive maximum Doppler shift and negative minimum Doppler shift, respectively. However, the reason why the energy is approximately uniform distributed in Case II is that the propagation scenario can be assumed to be a rich scattering environment in which the directions of incoming waves are uniform. In addition, the theoretical results fit well with the corresponding simulation ones obtained by (22).Figure 6
Theoretical TV-DPSDSpqt,f (Kpq = 0).Figure 7
Three typical Rx locations with respect to Tx.Figure 8
Theoretical and simulated TV-DPSDSpqt,f at three special cases (Kpq = 0).


### 3.4. TV-CCF

The absolute value of the TV-CCFρpq,p′q′t,δR can be computed using (16)-(19), as depicted in Figure 9. It can be easily observed that the TV-CCFs decrease as the antenna element spacing of the Rx increases. The TV-CCFs at different time instantst = 0 s, 1.5 s, 3 s (corresponding to xR = 0 m, 33.3 m, 66.6 m) are extracted in Figure 10, respectively. It is demonstrated that when the Rx gets closer to the Tx (fixed location, xT = 75 m), the TV-CCF decreases more sharply with the increasing of the antenna spacing. This trend is caused by the variation of multipath dispersion in the spatial domain. When the Rx is far from the Tx, the multipath is relatively concentrated in one direction, while it is getting more dispersed as the Rx closing to Tx. A lower CCF is friendly to the MIMO system working in the way of spatial diversity in which the sub-channels are expected to keep independent with each other. Moreover, the corresponding simulation results are computed utilizing (23) and the good correspondence between the theoretical and simulated results validates the correctness of the above derivations.Figure 9
Absolute values of the theoretical TV-CCFρpq,p′q′t,δR (Kpq = 0).Figure 10
Absolute values of the theoretical and simulated TV-CCFρpq,p′q′t,δR at different time instants (Kpq = 0).


## 3.1. Simulation Model

By assuming the infinite number of the scatterers (M,N,O→∞), the analytical model described in (2) is hard to realize. However, a corresponding simulation model with a finite number of scatterers can be developed. In terms of the model parameters, the discrete random values {xm,yn,zo} are generated by discrete uniform random parameters, while the other parameters are the same as the theoretical model, summarized in Table 2. To achieve a better approximation to the statistical properties of the theoretical model, the simulation results are obtained by averaging the results of multiple simulation runs (up to 104), and the number of scatterers per simulation run is selected as M×N×O=2000. Furthermore, the statistical properties of the simulation model are shown as (21)-(23). For TV-ACF, the simulation model is expressed as(21)r~pqSBt,τ=1Kpq+1·1MNO⁡∑m,n,o=1M,N,Oej2π∫t-τtfRmaxcos⁡αRm,nu-φvRcosβRm,n,ouduFor TV-DPSD, the simulation model is given as(22)S~pqSBt,f=∫-∞∞r~pqSBt,τe-j2πfτdτFor TV-CCF, i.e., the CCF at MSR, the simulation model is written as follows:(23)ρ~pq,p′q′SBt,δR=1Kpq+1·1MNO⁡∑m,n,o=1M,N,Oe-j2π/λδRq-q′cosϕRcosβRm,n,otcosγR-αRm,nt·e-j2π/λδRq-q′sinϕRsinβRm,n,otTable 2
Simulation parameters.
Parameters Value M T, MR 2, 2 f c (GHz) 1.8 v T, vR (km/h) 0, 80 (L,W,H) (m) (150, 5.2, 5) ( x T , y T , z T ) (m) (75, 2, 3) ( x R , y R , z R ) (m) (0, 0, 3) γ T, ϕT (rad) 0, 0 γ R, ϕR (rad) π / 2,  0 Ω p q (W) 1


## 3.2. TV-ACF

By using (12)-(14), the numerical results obtained for the TV-ACF rpqt,τ in NLOS scenario are shown in Figure 3. For clearly comparison purpose, the absolute values of the theoretical TV-ACF at different time instantst = 0 s, 1.5 s, 3 s (corresponding to xR = 0 m, 33.3 m, 66.6 m) are extracted as shown in Figure 4. Also, the corresponding simulation results can be obtained by (21). It can be seen that a good agreement between the theoretical and simulated results is achieved, which verifies the derivation of the proposed model. Furthermore, different variation trends of the ACF are resulted from the dynamic instantaneous phase which varies with the TV Doppler frequency. As the Rx gets closer to the Tx (fixed location, xT = 75 m), the channel autocorrelation behaves more sensitive to the time difference (see, for example, the curve att = 3 s). From Figure 5, it is obvious that the TV-ACF of the channel with the LOS component is higher than the case without the LOS component. In addition, the non-stationarity of channel characteristic is not notable in LOS dominated propagation condition.Figure 3
Absolute values of the theoretical TV-ACFrpqt,τ (Kpq = 0).Figure 4
Absolute values of the theoretical and simulated TV-ACFrpqt,τ at different time instants (Kpq = 0).Figure 5
Absolute values of the theoretical TV-ACFrpqt,τ with and without the LOS component at different time instants.


## 3.3. TV-DPSD

The theoretical TV-DPSD calculated by (15) is given in Figure 6. With the variation of the angle between the incoming wave and the direction of motion, the Doppler shift has a smooth drift from positive maximum to negative minimum. Figure 7 depicts three typical Rx locations with respect to Tx, the approach location (t = 0.45 s, xR = 10 m, Case I), the arrival location (t = 3.375 s, xR = 75 m, Case II), and the away location (t = 6.30 s, xR = 140 m, Case III). They are extracted from Figure 6 for further illustration as shown in Figure 8. A symmetric DPSD can be obtained in Case II while asymmetric DPSDs exist in Case I and Case III. In the non-isotropic environment like Case I and Case III, the AOA has a relatively concentrated distribution only in one direction. As a result, the energy focuses on the positive maximum Doppler shift and negative minimum Doppler shift, respectively. However, the reason why the energy is approximately uniform distributed in Case II is that the propagation scenario can be assumed to be a rich scattering environment in which the directions of incoming waves are uniform. In addition, the theoretical results fit well with the corresponding simulation ones obtained by (22).Figure 6
Theoretical TV-DPSDSpqt,f (Kpq = 0).Figure 7
Three typical Rx locations with respect to Tx.Figure 8
Theoretical and simulated TV-DPSDSpqt,f at three special cases (Kpq = 0).


## 3.4. TV-CCF

The absolute value of the TV-CCFρpq,p′q′t,δR can be computed using (16)-(19), as depicted in Figure 9. It can be easily observed that the TV-CCFs decrease as the antenna element spacing of the Rx increases. The TV-CCFs at different time instantst = 0 s, 1.5 s, 3 s (corresponding to xR = 0 m, 33.3 m, 66.6 m) are extracted in Figure 10, respectively. It is demonstrated that when the Rx gets closer to the Tx (fixed location, xT = 75 m), the TV-CCF decreases more sharply with the increasing of the antenna spacing. This trend is caused by the variation of multipath dispersion in the spatial domain. When the Rx is far from the Tx, the multipath is relatively concentrated in one direction, while it is getting more dispersed as the Rx closing to Tx. A lower CCF is friendly to the MIMO system working in the way of spatial diversity in which the sub-channels are expected to keep independent with each other. Moreover, the corresponding simulation results are computed utilizing (23) and the good correspondence between the theoretical and simulated results validates the correctness of the above derivations.Figure 9
Absolute values of the theoretical TV-CCFρpq,p′q′t,δR (Kpq = 0).Figure 10
Absolute values of the theoretical and simulated TV-CCFρpq,p′q′t,δR at different time instants (Kpq = 0).


## 4. Experimental Verification with Tunnel-Like Environment Measurement

A measurement campaign was carried out in a tunnel-like scenario at Zhongtian Technology Company, Nantong, China. As shown in Figure11, the overall length of the experimental tunnel is 100 m, with 50 m long rectangular part and 50 m long arched part. A concrete section exists in the junction which combines these two parts as shown in Figure 11(a). The 50 m long rectangular part has been considered for the measurement. The size of the rectangular cross section is 4.4 m (width)×3 m (height). The measurement was performed in three experiment regions denoted as A1, A2, and A3. In each area, a matrix of 1 × 8 receiving points was configured to simulate a virtual receiving array, as shown in Figure 11(b). The transmitting antenna was set at the middle of the tunnel, which was 0.5 m far away from the tunnel wall and 1.6 m high as illustrated in Figure 11(c). In the processing of the measured data, the channel impulse response (CIR) is directly generated by the cross-correlation of the measured data collected by the Rx with the replica of the transmitted PN sequence. Then the multipath components are detected with significant power contributions above the 5 dB of the estimated multipath noise floor. And finally, the narrowband CIR obtained by the complex sum over the delay domain is substituted into the expression (10) in [24] to extract the spatial correlation coefficient. As for the theoretical model, the spatial CCF can be obtained by (23) and the model parameters are set according to the measurement campaigns. Figure 12 shows that the theoretical curves can properly capture the measurement results for spatial CCF in A1, A2, and A3; especially a good performance is found in A2. This result demonstrates that the proposed model can generate a similar scattering environment as a realistic tunnel-like environment.Figure 11
Nantong short channel. (a) 3D tunnel geometric model. (b) Tunnel top view. (c) Tunnel sectional view.
(a) (b) (c)Figure 12
Comparison of the spatial CCF with various receive antenna spacings in (a) A1 area, (b) A2 area, and (c) A3 area.
(a) (b) (c)


## 5. Conclusion

In this paper, a 3D GBSM for non-stationary channels in tunnel environment has been proposed. By introducing the time-varying multipath length and AOA, the proposed model can be applied to simulate the non-stationary channels in real time. The different variation trends of TV-ACF at different time instants are related to the dynamic changes of the instantaneous phase in temporal. And it is obvious that the absolute values of TV-ACFs are higher in the propagation environment with the LOS components. Additionally, the TV-DPSDs for three typical cases have shown that the fluctuating trend of energy distribution in Doppler frequency arises from the variations of the scattering environment. It is also demonstrated that the TV-CCFs decrease with the increment in time, and when the multipath in the spatial domain is relatively dispersed, the decreasing rate is higher. Finally, a good fitting of the spatial CCFs of the theoretical model with the measured results verifies the feasibility of the proposed model.In future studies, the moving scatterers or the variation of the scatterers' distribution should be incorporated in the non-stationary model. For the dedicated cases in other tunnel scenarios, we will be devoted to establishing the applicable models to construct the spatial distribution of scatterers. Validation of the proposed channel model with other real-time measurement scenarios and with respect to other channel statistical properties is also an important topic for future research. Furthermore, the proposed channel model can also be studied under other signal frequency, tunnel curvature, and moving trajectory.



---

*Source: 6750153-2019-05-02.xml*