A 2D Non-Stationary Channel Model for Underwater Acoustic Communication Systems

Xiuming Zhu^1,2, Cheng-Xiang Wang^1,2,*, Ruofei Ma³
¹National Mobile Communications Research Laboratory, School of Information Science and Engineering,
Southeast University, Nanjing 210096, China.
²Purple Mountain Laboratories, Nanjing 211111, China.
³School of Information Science and Engineering, Harbin Institute of Technology, Weihai 264209, China.
^*Corresponding Author: Cheng-Xiang Wang
Email: {xm_zhu, chxwang}@seu.edu.cn, [email protected]

Abstract

Underwater acoustic (UWA) communication plays a key role in the process of exploring and studying the ocean. In this paper, a modified non-stationary wideband channel model for UWA communication in shallow water scenarios is proposed. In this geometry-based stochastic model (GBSM), multiple motion effects, time-varying angles, distances, clusters’ locations with the channel geometry, and the ultra-wideband property are considered, which makes the proposed model more realistic and capable of supporting long time/distance simulations. Some key statistical properties are investigated, including temporal autocorrelation function (ACF), power delay profile (PDP), average delay, and root mean square (RMS) delay spread. The impacts of multiple motion factors on temporal ACFs are analyzed. Simulation results show that the proposed model can mimic the non-stationarity of UWA channels. Finally, the proposed model is validated with measurement data.

Index Terms:

Underwater acoustic communication, shallow water, channel modeling, GBSM, non-stationarity

I Introduction

As a potential technology in supporting underwater communications in the sixth generation (6G) space-air-ground-sea integrated networks[1, 2], UWA communication plays a key role in exploring and studying the ocean. Since channel models are essential for the design and evaluation of communication systems[3], accurate UWA channel models with low complexity and good flexibility are indispensable.

Due to the complex propagation environments, UWA channels show several unique characteristics. UWA channels are affected by numerous motion factors[4, 5], which can be divided into: 1) intentional platform motion, e.g., autonomous underwater vehicle’s (AUV’s) vehicular motion; 2) unintentional drifting platform motion caused by the movement of water; 3) surface motion, which may fluctuate with waves and show cyclic patterns[6]. Because of the low speed of sound in the water (usually about 1500 m/s), the changes of transmission delays caused by motion effects cannot be neglected[7, 8]. Moreover, UWA channels may exhibit non-stationarity as a result of time-varying delays and ultra-wideband property, violating the wide sense stationary uncorrelated scattering (WSSUS) assumption[8, 9].

A number of studies have worked on the UWA communication channel modeling[10, 11, 4, 12, 5, 13, 14, 15]. However, no standardized channel model has been proposed yet[16]. In [10, 11], typical solutions were used to calculate the acoustic field, e.g., ray tracing and parabolic approximation. However, these models, while accurate, are deterministic and lack of flexibility. In [4, 12, 5, 13, 14, 15], several stochastic channel models were proposed to characterize UWA channels. In these models, GBSMs have good balance among accuracy, complexity, and flexibility[13, 14, 15], and are widely used in the modeling of wireless communication channels. In [13, 14], GBSMs based on WSSUS assumption were proposed and analyzed, which cannot capture the non-stationarity of UWA channels. In [15], the effects of platform and scatterers’ motion were considered. However, the velocities of the transmitter (Tx) and receiver (Rx) were assumed to be constant, which could not model the random drifting motion in long time simulations. The delay changes caused by surface motion were modeled as Gaussian processes, which cannot characterize the cyclic patterns in channel properties. Moreover, the model neglected the change of the channel geometry due to significant displacements of Tx and Rx, thus making the model not suitable for long time/distance simulations, e.g., AUV’s application scenarios.

To the best of the authors’ knowledge, non-stationary GBSMs for shallow water scenarios considering multiple motion factors and allowing for long time/distance communication scenarios are still missing in the literature. This work is aiming to fill the research gap. The major contributions and novelties of this work can be summarized as follows.

1.

The proposed model considers two kinds of platform motion (constant intentional motion and random drifting motion) and sinusoidal motion of surface scatterers. The distances, angles, especially the clusters’ locations are modeled as time-varying parameters.
2.

The time-frequency varying transmisson losses caused by time-varying propagation distances and underwater frequency-dependent absorption loss are considered in the proposed model.
3.

The proposed twin-cluster GBSM supports single-bounce (SB) and multiple-bounce (MB) propagation. For MB propagation, the angles-of-arrival (AoAs) are independent with the angles-of-departure (AoDs), while these angles are correlated with geometric relationships for SB propagation.
4.

Based on the proposed model, some important statistical properties such as ACF, PDP, average delay, and RMS delay spread are studied and analyzed. The proposed GBSM is also validated with measurement data.

The rest of this paper is organized as follows. In Section II, a modified non-stationary shallow water GBSM is introduced in detail. The key statistical properties of the model are derived in Section III. Section IV presents the simulation results and analysis. Conclusions are finally drawn in Section V.

II A Modified Non-stationary GBSM for Shallow Water UWA Communication

Refer to caption — Figure 1: Illustration of different propagation components ( $N_{S}=1,N_{B}=1$ ).

Let us consider a single-input single-output (SISO) communication system in a shallow water scenario, where both Tx and Rx are located in the water. The sound speed in the water is related to salinity, temperature, and pressure. In shallow water scenarios, these factors are usually constant because of the mixing of wind, so is the sound speed ( $\approx 1500$ m/s)[7, 17]. Thus, acoustic signals can be assumed to propagate along straight lines in shallow water scenarios. The acoustic signals may reach the Rx directly or reflect several times between two boundaries. Fig. 1 gives an illustration of different propagation components, i.e., line-of-sight (LoS), downward arrival (DA), and upward arrival (UA) paths. Due to the roughness of the boundaries, it is assumed that each DA (UA) path is comprised of $N_{s\tilde{b}}^{\rm DA}$ ( $M_{b\tilde{s}}^{\rm UA}$ ) diffuse scattering rays. The specular reflection ray is called the macro-scattering ray, while the diffuse rays are defined as micro-scattering rays [13]. The micro-scatterers are assumed to be clustered around the macro-cluster, i.e., the specular reflection point. Each macro-scattering ray is assumed to be the average of micro-scattering rays. Assuming that each DA (UA) path will contact the surface (bottom) at most $N_{S}$ ( $N_{B}$ ) time(s) with $s$ ( $\tilde{s}$ ) surface interaction(s) and $\tilde{b}$ ( $b$ ) bottom interaction(s), then there are constraints $1\leq s\leq N_{S}$ ( $1\leq b\leq N_{B}$ ) and $s-1\leq\tilde{b}\leq s$ ( $b-1\leq\tilde{s}\leq b$ ).

Fig. 2 shows the proposed modified GBSM. For clarity, only the LoS path and the $m$ th ray in the UA path ( $b=1$ , $\tilde{s}=1$ ) are illustrated. Considering a MB propagation, the AoDs of rays are only related to the first-bounce cluster $C_{b\tilde{s}}^{{\rm UA},A}$ while the AoAs of rays are only related to the last-bounce cluster $C_{b\tilde{s}}^{{\rm UA},Z}$ [18]. When $b+\tilde{s}=1$ is satisfied, the MB propagation is reduced to SB propagation. The intentional motion velocity of Tx (Rx) is assumed to be constant, with speed $V_{M}^{T}$ ( $V_{M}^{R}$ ) and travel angle $\alpha_{M}^{T}$ ( $\alpha_{M}^{R}$ ). The random drifting movement of Tx (Rx) is described by the displacement $d_{D}^{T}(t)$ ( $d_{D}^{R}(t)$ ) and travel angle $\alpha_{D}^{T}(t)$ ( $\alpha_{D}^{R}(t)$ ). The $m$ th scatterer on the surface has speed $V_{b\tilde{s},m}^{S,{\rm UA},A}(t)$ ( $V_{b\tilde{s},m}^{S,{\rm UA},Z}(t)$ ) with constant travel angle $\alpha^{S}$ . The water depth, Tx depth, Rx depth, and the horizontal distance between Tx and Rx are denoted by $h_{S}$ , $h_{T}(t)$ , $h_{R}(t)$ , and $D(t)$ , respectively. The transmission distance, AoD, and AoA of the LoS path are denoted by $d^{\rm LoS}(t)$ , $\phi_{\rm LoS}^{T}(t)$ and $\phi_{\rm LoS}^{R}(t)$ , respectively. The AoD and AoA of the $m$ th ray in the UA path are denoted by $\phi_{b\tilde{s},m}^{{\rm UA},T}(t)$ and $\phi_{b\tilde{s},m}^{{\rm UA},R}(t)$ , respectively. The transmission distances Tx $-S_{b\tilde{s},m}^{{\rm UA},A}$ , $S_{b\tilde{s},m}^{{\rm UA},A}-S_{b\tilde{s},m}^{{\rm UA},Z}$ , and $S_{b\tilde{s},m}^{{\rm UA},Z}-$ Rx are denoted by $d_{b\tilde{s},m}^{{\rm UA},T}(t)$ , $d_{b\tilde{s},m}^{{\rm UA},S}(t)$ , and $d_{b\tilde{s},m}^{{\rm UA},R}(t)$ , respectively. For DA paths, there are similar denotations.

II-A Channel Transfer Function (CTF)

Based on the proposed GBSM shown in Fig. 2, the CTF of the channel model is the superposition of LoS, DA, and UA components, which can be expressed as

	$\displaystyle H(t,f)$	$\displaystyle=\sqrt{\frac{K}{K+1}}H^{\rm LoS}(t,f)$		(1)
		$\displaystyle+\sqrt{\frac{\eta_{\rm DA}}{K+1}}H^{\rm DA}(t,f)+\sqrt{\frac{\eta_{\rm UA}}{K+1}}H^{\rm UA}(t,f)$		(1)

where $K$ is the Rice factor, $\eta_{\rm DA}$ ( $\eta_{\rm UA}$ ) denotes the ratio of DA (UA) component’s power to the total non-LoS (NLoS) power, with constraint $\eta_{\rm DA}+\eta_{\rm UA}=1$ . The three propagation components in (1) are modeled as

$H^{\rm LoS}(t,f)=a^{\rm LoS}(t,f)e^{-j2\pi\left(f_{c}+f\right)\tau^{\rm LoS}(t)}$		(2a)

$\displaystyle H^{\rm DA}(t,f)=$	$\displaystyle\frac{1}{\sqrt{2N_{S}N_{s\tilde{b}}^{\rm DA}}}\sum\nolimits_{s=1}^{N_{S}}\sum\nolimits_{\tilde{b}=s-1}^{s}\sum\nolimits_{n=1}^{N_{s\tilde{b}}^{\rm DA}}$	(2b)
	$\displaystyle a_{s\tilde{b},n}^{\rm DA}(t,f)e^{j\theta_{s\tilde{b},n}^{\rm DA}-j2\pi\left(f_{c}+f\right)\tau_{s\tilde{b},n}^{\rm DA}(t)}$	(2b)

$\displaystyle H^{\rm UA}(t,f)=$	$\displaystyle\frac{1}{\sqrt{2N_{B}M_{b\tilde{s}}^{\rm UA}}}\sum\nolimits_{b=1}^{N_{B}}\sum\nolimits_{\tilde{s}=b-1}^{b}\sum\nolimits_{m=1}^{M_{b\tilde{s}}^{\rm UA}}$	(2c)
	$\displaystyle a_{b\tilde{s},m}^{\rm UA}(t,f)e^{j\theta_{b\tilde{s},m}^{\rm UA}-j2\pi\left(f_{c}+f\right)\tau_{b\tilde{s},m}^{\rm UA}(t)}.$	(2c)

Here, $f_{c}$ is the carrier frequency, $a^{\rm LoS}(t,f)$ and $\tau^{\rm LoS}(t)$ are path gain and transmission delay of the LoS path, respectively. Parameters $a_{s\tilde{b},n}^{\rm DA}(t,f)$ ( $a_{b\tilde{s},m}^{\rm UA}(t,f)$ ), $\tau_{s\tilde{b},n}^{\rm DA}(t)$ ( $\tau_{b\tilde{s},m}^{\rm UA}(t)$ ), and $\theta_{s\tilde{b},n}^{\rm DA}$ ( $\theta_{b\tilde{s},m}^{\rm UA}$ ) are the gain, delay, and initial phase shift of the $n$ th ( $m$ th) ray in the DA (UA) path, respectively. The initial phase shifts are modeled as independent random variables with uniform distributions over $[0,2\pi)$ . Derivations of the rest model parameters are given in the remainder of this section.

II-B Time-varying Delays

Due to multiple motion factors in UWA channels, transmission delays are changing with time and the motion factors can be divided into three types[4, 5]: 1) intentional platform motion, e.g., AUV’s vehicular motion; 2) unintentional drifting platform motion; 3) surface motion. The speeds $V_{M}^{X}$ and directions of travel $\alpha_{M}^{X}$ of intentional motion are often given in measurement campaigns and are assumed to be constant here. For the drifting motion, considering the randomness of the water movement, it is assumed that the drifting velocity vectors $\vec{V}_{D}^{X}(t)$ remain unchanged for a period of time, and change at intervals with the speeds randomly generated between $[v_{D}^{\rm min},v_{D}^{\rm max}]$ and the angles of travel randomly generated between $[0,2\pi)$ . The change frequency is denoted by $f_{v}^{D}$ . Then the lengths of displacements and directions of drifting are calculated as $d_{D}^{X}(t)=\|\int\vec{V}_{D}^{X}(t)dt\|$ and $\alpha_{D}^{X}(t)=\arg\{\int\vec{V}_{D}^{X}(t)dt\}$ , respectively. Here, $X=\{T,R\}$ , where the superscripts $T$ and $R$ denote the Tx and Rx, respectively. For the surface motion, the relative motion speed of the $k$ th scatterer on the surface to the corresponding cluster is modeled as $V_{k}^{S}(t)=2\pi f^{S}A^{S}{\cos}(2\pi f^{S}t+\theta_{k}^{S})$ , with a constant angle of travel $\alpha^{S}$ . Similar assumption can be found in [4]. The phase term $\theta_{k}^{S}$ is randomly generated between $[0,2\pi)$ , characterizing the randomness of the initial periodic motion of each scatterer on the surface.

II-B1 LoS

For the LoS path, the time-varying delay is calculated as $\tau^{\rm LoS}(t)=d^{\rm LoS}(t)/c$ , where $c$ denotes the sound speed in the water. The time-varying distance is given by

$\displaystyle d^{\rm LoS}(t)\approx$	$\displaystyle\sqrt{D^{2}(t)+\left[h_{R}(t)-h_{T}(t)\right]^{2}}$	(3)
	$\displaystyle-d_{D}^{T}(t)\cos\left[\alpha_{D}^{T}(t)-\phi_{\rm LoS}^{T}(t)\right]$
	$\displaystyle-d_{D}^{R}(t)\cos\left[\phi_{\rm LoS}^{R}(t)-\alpha_{D}^{R}(t)\right].$

Note that the impacts of drifting motion on the positions of Tx and Rx are small, so only the influence of intentional motion on the channel geometry is considered. Denoting the initial distances as $D(t_{0})$ , $h_{T}(t_{0})$ , and $h_{R}(t_{0})$ , the time-varying distances in (3) are decided by

	$D(t)=D\left(t_{0}\right)-V_{M}^{T}t\cdot\cos\alpha_{M}^{T}+V_{M}^{R}t\cdot\cos\alpha_{M}^{R}$		(4a)
	$h_{T}(t)=h_{T}\left(t_{0}\right)+V_{M}^{T}t\cdot\sin\alpha_{M}^{T}$		(4b)
	$h_{R}(t)=h_{R}\left(t_{0}\right)+V_{M}^{R}t\cdot\sin\alpha_{M}^{R}.$		(4c)

The AoD and AoA of LoS path in (3) can be calculated as $\phi_{\rm Los}^{T}(t)=\arctan\frac{h_{R}(t)-h_{T}(t)}{D(t)}$ and $\phi_{\rm LoS}^{R}(t)=\phi_{\rm LoS}^{T}(t)+\pi$ , respectively.

II-B2 NLoS (DA & UA)

Based on geometric relationships, the transmission distances and characteristic angles of incidence (AOIs) of macro-scattering rays are determined by

$\displaystyle d_{s\tilde{b}}^{\rm DA}(t)$	$\displaystyle=\sqrt{D^{2}(t)+\left[2sh_{S}+(-1)^{s-\tilde{b}}h_{T}(t)-h_{R}(t)\right]^{2}}$	(5)
$\displaystyle d_{b\tilde{s}}^{\rm UA}(t)$	$\displaystyle=\sqrt{D^{2}(t)+\left[2\tilde{s}h_{S}-(-1)^{b-\tilde{s}}h_{T}(t)+h_{R}(t)\right]^{2}}$	(6)
$\displaystyle\varphi_{s\tilde{b}}^{\rm DA}(t)$	$\displaystyle=\arctan\frac{D(t)}{2sh_{S}+(-1)^{s-\tilde{b}}h_{T}(t)-h_{R}(t)}$	(7)
$\displaystyle\varphi_{b\tilde{s}}^{\rm UA}(t)$	$\displaystyle=\arctan\frac{D(t)}{2\tilde{s}h_{S}-(-1)^{b-\tilde{s}}h_{T}(t)+h_{R}(t)}$	(8)

where examples of AOIs $\varphi_{s\tilde{b}}^{\rm DA}$ and $\varphi_{b\tilde{s}}^{\rm UA}$ are shown in Fig. 1.

For brevity, we use simplified symbols to describe the derivations. Considering the propagation route Tx $-S^{A}-S^{Z}-$ Rx of a ray in a NLoS path, where $S^{A}$ ( $S^{Z}$ ) denotes the scatterer in the first- (last-) bounce cluster $C^{A}$ ( $C^{Z}$ ), the transmission delay of the ray can be obtained by $\tau_{r}(t)=[d_{r}^{T}(t)+d_{r}^{S}(t)+d_{r}^{R}(t)]/c$ , where $d_{r}^{T}(t)$ , $d_{r}^{S}(t)$ , and $d_{r}^{R}(t)$ denote $d({\rm Tx}-S^{A})$ , $d(S^{A}-S^{Z})$ , and $d(S^{Z}-{\rm Rx})$ , respectively. The AoD and AoA of the micro-scattering ray are denoted by $\phi_{r}^{T}(t)$ and $\phi_{r}^{R}(t)$ , respectively. Similarly, $d_{c}^{T}(t)$ , $d_{c}^{S}(t)$ , and $d_{c}^{R}(t)$ denote $d({\rm Tx}-C^{A})$ , $d(C^{A}-C^{Z})$ , and $d(C^{Z}-{\rm Rx})$ , respectively. The transmission distance, AoD, AoA, and characteristic AOI of the macro-scattering ray are denoted by $d_{c}(t)$ , $\phi_{c}^{T}(t)$ , $\phi_{c}^{R}(t)$ , and $\varphi_{c}(t)$ , respectively. Note that the relevant distances need to be classified and calculated separately according to different interactions of rays with the surface and the bottom. For DA (UA) paths, $C^{Z}$ s are located on the surface (bottom), $C^{A}$ s are on the bottom (surface) if $s=\tilde{b}$ ( $b=\tilde{s}$ ) and on the surface (bottom) if not.

For MB propagation (DA paths: $s+\tilde{b}\neq 1$ , UA paths: $b+\tilde{s}\neq 1$ ), the first and last distance can be derived as

	$d_{r}^{T}(t)\approx\left\{\begin{array}[]{l}A_{r}^{T}(t)+\frac{h_{S}-h_{T}(t)}{\sin\phi_{r}^{T}(t)}-B_{r}^{T}(t),C^{A}\text{ at surface }\\ \frac{h_{T}(t)}{\sin\left[2\pi-\phi_{r}^{T}(t)\right]}-B_{r}^{T}(t),C^{A}\text{ at bottom }\end{array}\right.$		(9a)
	$d_{r}^{R}(t)\approx\left\{\begin{array}[]{l}A_{r}^{R}(t)+\frac{h_{S}-h_{R}(t)}{\sin\left[\pi-\phi_{r}^{R}(t)\right]}-B_{r}^{R}(t),C^{Z}\text{ at surface }\\ \frac{h_{R}(t)}{\sin\left[\phi_{r}^{R}(t)-\pi\right]}-B_{r}^{R}(t),C^{Z}\text{ at bottom}.\end{array}\right.$		(9b)

In (9), $A_{r}^{X}(t)=A^{S}\sin\left(2\pi f^{S}t+\theta_{r}^{S,X}\right)\cos\left[\phi_{r}^{X}(t)-\alpha^{S}\right]$ , $B_{r}^{X}(t)=d_{D}^{X}(t)\cos\left[\alpha_{D}^{X}(t)-\phi_{r}^{X}(t)\right]$ , $\phi_{r}^{X}(t)$ are assumed to have Gaussian distributions, and the mean AoA and AoD are determined by $\phi_{c}^{X}(t)$ , where $X=\{T,R\}$ . The $\theta_{r}^{S,T}$ ( $\theta_{r}^{S,R}$ ) denotes the phase of sinusoidal speed of the surface scatterer in $C^{A}$ ( $C^{Z}$ ). Two parameters are used to control the angle spreads, i.e., $\phi_{r}^{X}(t)\sim N(\phi_{c}^{X}(t),\sigma_{\phi,s}^{2})$ for surface-interacting rays and $\phi_{r}^{X}(t)\sim N(\phi_{c}^{X}(t),\sigma_{\phi,b}^{2}$ ) for bottom-interacting rays. For $C^{A}$ ( $C^{Z}$ ), $\phi_{c}^{T}(t)=\pi/2-\varphi_{c}(t)$ ( $\phi_{c}^{R}(t)=\pi/2+\varphi_{c}(t)$ ) and $\phi_{c}^{T}(t)=3\pi/2+\varphi_{c}(t)$ ( $\phi_{c}^{R}(t)=3\pi/2-\varphi_{c}(t)$ ) when interacting with surface and bottom, respectively.

The second distance of the micro-scattering ray is assumed to be very close to the distance of the macro-scattering ray, and is modeled as $d_{r}^{S}(t)=d_{c}^{S}(t)\cdot e^{\Delta d^{S}}$ , where the $\Delta d^{S}$ is a zero-mean Gaussian random variable with variance $\sigma_{ds}^{2}$ . Parameter $\sigma_{ds}$ is used to characterize small differences between the second distances of micro-scattering rays and the macro-scattering ray and it is set as a small value of 0.001 arbitrarily in this paper. Then there is $d_{c}^{S}(t)=d_{c}(t)-d_{c}^{T}(t)-d_{c}^{R}(t)$ , where $d_{c}^{X}(t)=H(t)/\cos(\varphi_{c}(t)$ , $X=\{T,R\}$ . For surface-interacting (bottom-interacting) clusters, $H(t)=h_{S}-h_{X}(t)$ ( $H(t)=h_{X}(t)$ ), where $X=T$ for $C^{A}$ and $X=R$ for $C^{Z}$ .

For SB propagation (DA path: $s=1,\tilde{b}=0$ , UA path: $b=1,\tilde{s}=0$ ) the angles are correlated in terms of geometric relationships and the second distance is zero, i.e., $d_{r}^{S}(t)=0$ . The AoDs of DA and UA paths in SB propagation can be determined by $\phi_{r}^{T}(t)=\arctan\frac{h_{S}-h_{T}(t)}{D(t)-\left[h_{S}-h_{R}(t)\right]/\tan\left[\pi-\phi_{r}^{R}(t)\right]}$ and $\phi_{r}^{T}(t)=2\pi-\arctan\frac{h_{T}(t)}{D(t)-h_{R}(t)/\tan\left[\phi_{r}^{R}(t)-\pi\right]}$ , respectively, where $\theta_{r}^{S,T}=\theta_{r}^{S,R}$ for DA path.

II-C Time-frequency Varying Gains

II-C1 LoS

The gain of the LoS path can be expressed as

a^{\rm LoS}(t,f)=L_{s}\left(d^{\rm LoS}(t)\right)L_{a}\left(d^{\rm LoS}(t),f\right)

(10)

where the geometric spreading loss coefficient at distance $d$ (m) is $L_{s}(d)=1/d$ [15] as we assume that Tx is a point source approximately which generates spherical spreading. The absorption loss coefficient in (10) can be expressed as $L_{a}(d,f)=10^{-\frac{d\cdot\alpha(f)}{20000}}$ [15], where $f$ is the frequency (kHz), $\alpha(f)$ is the frequency-dependent absorption parameter (dB/km) which can be given by the empirical Thorp model $\alpha(f)=\frac{0.11f^{2}}{1+f^{2}}+\frac{44f^{2}}{4100+f^{2}}+2.75\times 10^{-4}f^{2}+0.003$ [19].

II-C2 NLoS (DA & UA)

The gains of DA (UA) paths are related to the transmission distances, AOIs and frequency. As the distances and angles of micro-scattering rays are very close to that of the corresponding macro-scattering ray, we assume the gains of micro-scattering rays in a NLoS path are equal to that of the macro-scattering ray, i.e., $a_{s\tilde{b},n}^{\rm DA}(t,f)=a_{s\tilde{b}}^{\rm DA}(t,f)$ and $a_{b\tilde{s},m}^{\rm UA}(t,f)=a_{b\tilde{s}}^{\rm UA}(t,f)$ for the sake of simplicity. The gains can be expressed as

	$\displaystyle a_{s\tilde{b}}^{\rm DA}(t,f)$	$\displaystyle=L_{s}\left(d_{s\tilde{b}}^{\rm DA}(t)\right)L_{a}\left(d_{s\tilde{b}}^{\rm DA}(t),f\right)L_{b}\left(\varphi_{s\tilde{b}}^{\rm DA}(t)\right)^{\tilde{b}}$		(11)
	$\displaystyle a_{b\tilde{s}}^{\rm UA}(t,f)$	$\displaystyle=L_{s}\left(d_{b\tilde{s}}^{\rm UA}(t)\right)L_{a}\left(d_{b\tilde{s}}^{\rm UA}(t),f\right)L_{b}\left(\varphi_{b\tilde{s}}^{\rm UA}(t)\right)^{b}$		(12)

where $L_{b}(\cdot)$ is the bottom reflection loss and is given by [19]

L_{b}(\varphi)=\left|\frac{\left(\rho_{b}/\rho_{w}\right)\cos(\varphi)-\sqrt{\left(c_{w}/c_{b}\right)^{2}-\sin^{2}(\varphi)}}{\left(\rho_{b}/\rho_{w}\right)\cos(\varphi)+\sqrt{\left(c_{w}/c_{b}\right)^{2}-\sin^{2}(\varphi)}}\right|

(13)

where $\varphi$ is the AoI at the bottom, $\rho_{b}$ ( $\rho_{w}$ ) denotes the density of the bottom (water), and $c_{b}$ ( $c_{w}$ ) denotes the sound speed in the bottom (water).

III Statistical properties

III-A Time-Frequency Correlation Function (TF-CF)

The TF-CF of the channel model can be defined as

R_{H}(t,f;\Delta t,\Delta f)=\mathrm{E}\left\{H(t,f)H^{*}(t-\Delta t,f-\Delta f)\right\}

(14)

where $\mathrm{E}\{\cdot\}$ denotes the ensemble average and $(\cdot)^{*}$ is the complex conjugate operation. Assuming that LoS, DA, and UA paths are independent with each other, by substituting (1) and (2) into (14), the TF-CF can be expressed as

\begin{array}[]{l}R_{H}(t,f;\Delta t,\Delta f)=\frac{K}{K+1}R_{H}^{\rm LoS}(t,f;\Delta t,\Delta f)\\ \quad+\frac{\eta_{\rm DA}}{2N_{S}(K+1)}\sum_{s=1}^{N_{S}}\sum_{\tilde{b}=s-1}^{s}R_{H,s\tilde{b}}^{\rm DA}(t,f;\Delta t,\Delta f)\\ \quad+\frac{\eta_{\rm UA}}{2N_{B}(K+1)}\sum_{b=1}^{N_{B}}\sum_{\tilde{s}=b-1}^{b}R_{H,b\tilde{s}}^{\rm UA}(t,f;\Delta t,\Delta f).\end{array}

(15)


$\displaystyle R_{H}^{\rm LoS}(t,f;\Delta t,\Delta f)$	$\displaystyle=\mathrm{E}\left\{A^{\rm LoS}e^{-j2\pi\left(f_{c}+f\right)\left[\tau^{\rm LoS}(t)-\tau^{\rm LoS}(t-\Delta t)\right]-j2\pi\Delta f\tau^{\rm LoS}(t-\Delta t)}\right\}$	(16a)
$\displaystyle R_{H,s\tilde{b}}^{\rm DA}(t,f;\Delta t,\Delta f)$	$\displaystyle=\mathrm{E}\left\{\frac{1}{N_{s\tilde{b}}^{\rm DA}}\sum_{n=1}^{N_{s\tilde{b}}^{\rm DA}}A_{s\tilde{b},n}^{\rm DA}e^{-j2\pi\left(f_{c}+f\right)\left[\tau_{s\tilde{b},n}^{\rm DA}(t)-\tau_{s\tilde{b},n}^{\rm DA}(t-\Delta t)\right]-j2\pi\Delta f\tau_{s\tilde{b},n}^{\rm DA}(t-\Delta t)}\right\}$	(16b)
$\displaystyle R_{H,b\tilde{s}}^{\rm UA}(t,f;\Delta t,\Delta f)$	$\displaystyle=\mathrm{E}\left\{\frac{1}{M_{b\tilde{s}}^{\rm UA}}\sum_{m=1}^{M_{b\tilde{s}}^{\rm UA}}A_{b\tilde{s},m}^{\rm UA}e^{-j2\pi\left(f_{c}+f\right)\left[\tau_{b\tilde{s},m}^{\rm UA}(t)-\tau_{b\tilde{s},m}^{\rm UA}(t-\Delta t)\right]-j2\pi\Delta f\tau_{b\tilde{s},m}^{\rm UA}(t-\Delta t)}\right\}$	(16c)

Denoting that $A^{\rm LoS}=a^{\rm LoS}(t,f)a^{\rm LoS}(t-\Delta t,f-\Delta f)$ , $A_{s\tilde{b},n}^{\rm DA}=a_{s\tilde{b},n}^{\rm DA}(t,f)a_{s\tilde{b},n}^{\rm DA}(t-\Delta t,f-\Delta f)$ , and $A_{b\tilde{s},m}^{\rm UA}=a_{b\tilde{s},m}^{\rm UA}(t,f)a_{b\tilde{s},m}^{\rm UA}(t-\Delta t,f-\Delta f)$ , the correlation functions of LoS, DA, and UA components are given by (16a) $-$ (16c), shown at the top of the next page. The temporal ACF can be obtained by the TF-CF, i.e., $R_{H}^{\rm ACF}(t,f;\Delta t)=R_{H}(t,f;\Delta t,0)$ .

III-B PDP

The time-frequency varying PDP can be expressed as

\begin{array}[]{l}P(t,f;\tau)=\frac{K}{K+1}P_{\rm LoS}(t,f;\tau)\\ +\frac{\eta_{\rm DA}}{2N_{S}N_{s\tilde{b}}^{\rm DA}(K+1)}\sum_{s=1}^{N_{S}}\sum_{\tilde{b}=s-1}^{s}\sum_{n=1}^{N_{s\tilde{b}}^{\rm DA}}P_{s\tilde{b},n}^{\rm DA}(t,f;\tau)\\ +\frac{\eta_{\rm UA}}{2N_{B}M_{b\tilde{s}}^{\rm UA}(K+1)}\sum_{b=1}^{N_{B}}\sum_{\tilde{s}=b-1}^{b}\sum_{m=1}^{M_{b\tilde{s}}^{\rm UA}}P_{b\tilde{s},m}^{\rm UA}(t,f;\tau).\end{array}

(17)

In (17), the PDPs of three propagation components are determined as $P_{\rm LoS}(t,f;\tau)=\left[a_{\rm LoS}(t,f)\right]^{2}\delta\left(\tau-\tau_{\rm LoS}(t)\right)$ , $P_{s\tilde{b},n}^{\rm DA}(t,f;\tau)=\left[a_{s\tilde{b},n}^{\rm DA}(t,f)\right]^{2}\delta\left(\tau-\tau_{s\tilde{b},n}^{\rm DA}(t)\right)$ , and $P_{b\tilde{s},m}^{\rm UA}(t,f;\tau)=\left[a_{b\tilde{s},m}^{\rm UA}(t,f)\right]^{2}\delta\left(\tau-\tau_{b\tilde{s},m}^{\rm UA}(t)\right)$ , respectively.

III-C Average Delay and RMS Delay Spread

The average delay and the RMS delay spread can be obtained by the PDP, and can be calculated as

\mu_{\tau}(t,f)=\frac{\sum_{\tau}\tau\cdot P(t,f;\tau)}{\sum_{\tau}P(t,f;\tau)}

(18)

\sigma_{\tau}(t,f)=\sqrt{\frac{\sum_{\tau}\left(\tau-\mu_{\tau}(t,f)\right)^{2}\cdot P(t,f;\tau)}{\sum_{\tau}P(t,f;\tau)}}.

(19)

IV Results and Analysis

In this section, some key statistical properties of the proposed model are simulated and verified with the corresponding measurement data. In the simulation, we consider a scenario with geometry parameters given as $D(t_{0})=2000$ m, $h_{S}=100$ m, $h_{T}(t_{0})=50$ m, and $h_{R}(t_{0})=80$ m. The sound speed $c(c_{w})$ is set as $1500$ m/s [7]. We assume the drifting motion changes once per second and consider vertical surface motion, i.e., $f_{v}^{D}=1$ Hz and $\alpha^{S}=\pi/2$ . Since shallow UWA communication usually accounts for long-range communication, the transmission distance is usually large compared with the water depth. Thus, we assume that the angle spread is small and set the corresponding parameter as $\sigma_{\phi,s}=\sigma_{\phi,b}=0.015$ . And refer to [13], we set $N_{B}=N_{S}=2$ , $\eta_{\rm DA}=\eta_{\rm UA}=0.5$ , $\rho_{b}/\rho_{w}=1.5$ , and $c_{b}=1600$ m/s.

Fig. 3 illustrates the normalized absolute values of temporal ACFs of the proposed model under different motion influences with and without LoS component. It can be observed that multiple motion factors in UWA channels make the temporal ACFs decay faster. The temporal correlation remains higher and is less affected by motion factors with the LoS component. It can also be noticed that as the amplitude of surface scatterers’ periodic motion ( $A^{S}$ ) increases, the ACFs decay faster. Note that the surface motion speed is proportional to $A^{S}$ , we can conclude that larger motion speed will attenuate the temporal correlation more dramatically. It can also be predicted that UWA channel coherence time will be lower when ships pass through or the weather is windy. Moreover, the analytical results have a good consistency with the simulated results, validating the correctness of the derivations and simulations.

Fig. 4 presents the temporal ACFs at different time instants and different carrier frequencies. We can observe different temporal ACFs at $t=0$ s, $5$ s, and $10$ s, which are resulted from time-varying angles, distances, and clusters’ locations in the proposed model, illustrating that our model can mimic the non-stationarity of UWA channels. Besides, the distinction between temporal ACFs at different carrier frequencies shows that the magnitude of the coherence time is related to the communication frequency band. The higher the carrier frequency, the shorter the coherence time. In addition, the simulated results of ACFs can fit the analytical results well.

Fig. 5 shows the PDPs of the proposed model by using the mean gains and delays of rays within clusters. All transmission delays are normalized with respect to the delay of the first arrival ray. We can see that the PDPs are different at $t=0$ s and $5$ s, and have differences at $f_{c}=15$ kHz and $100$ kHz, which shows the non-stationarity of the proposed model again. The strongest ray is the first arrival ray in the isovelocity shallow water scenarios as can be seen in Fig. 5.

TABLE I: Comparison of Delay Statistics
(

K=1.44

\sigma_{\phi,s}=\sigma_{\phi,b}=0.02317

\eta_{\rm DA}=\eta_{\rm UA}=0.5

N_{S}=N_{B}=1

Delay Statistics	Measurement data in [13]	Simulation model in [13]	The proposed model
Average Delay	$1.5$ ms	$1.491$ ms	$1.505$ ms
RMS Delay Spread	$2.4$ ms	$2.410$ ms	$2.399$ ms

Table I gives the comparison of delay statistics of the simulation model in [13] and the proposed model with the corresponding measurement data in [13]. The parameters of measurement campaign are given in as $D(t_{0})=1500$ m, $h_{S}=80$ m, $h_{T}(t_{0})\approx 34.5$ m, $h_{R}(t_{0})\approx 36$ m, $c=1440$ m/s, $V_{M}^{T}=V_{M}^{R}=v_{D}^{\rm min}=v_{D}^{\rm max}=0$ m/s, $f_{c}=17$ kHz, $c_{b}/c_{w}=1.11$ , and $\rho_{b}/\rho_{w}=1.5$ [13] . The weather of the measurement scenario in [13] was rainy and windy, so we set parameters of surface motion as $A^{S}=2$ , $f^{S}=0.1$ Hz, and $\alpha^{S}=\pi/2$ [4]. With the rest of parameters chosen according to the estimation procedure introduced in [20] which is based on the minimum mean square error criterion, the proposed model matches better with measurement data, showing the practicality of our model.

V Conclusions

In this paper, a modified non-stationary GBSM for shallow water UWA communication systems has been proposed. Some important statistical properties have been investigated, and the analytical results have a good consistency with the simulated results, validating the correctness of the derivations and simulations. Simulation results have illustrated that multiple motion factors have great influence on temporal ACFs, so as on UWA communication systems. The impact of unintentional motion factors cannot be ignored. The fact that simulation results of the proposed model vary with time and frequency has shown that our model can mimic the non-stationarity of UWA channels. In addition, simulation results of statistical properties fit well with measurements, illustrating the usefulness of the proposed model. In our future work, we will try to consider the change of sound speed to improve the applicability of the channel model in deep waters.

Acknowledgment

This work was supported by the National Key R&D Program of China under Grant 2018YFB1801101, the National Natural Science Foundation of China (NSFC) under Grant 61960206006, the Frontiers Science Center for Mobile Information Communication and Security, the High Level Innovation and Entrepreneurial Research Team Program in Jiangsu, the High Level Innovation and Entrepreneurial Talent Introduction Program in Jiangsu, the Research Fund of National Mobile Communications Research Laboratory, Southeast University, under Grant 2020B01, the Fundamental Research Funds for the Central Universities under Grant 2242020R30001, and the EU H2020 RISE TESTBED2 project under Grant 872172.

References

[1] C.-X. Wang, J. Huang, H. Wang, X. Gao, X.-H. You, and Y. Hao, “6G wireless channel measurements and models: Trends and challenges,” IEEE Veh. Technol. Mag., vol. 15, no. 4, pp. 22–32, Dec. 2020.
[2] X.-H. You, C.-X. Wang, J. Huang, et al., “Towards 6G wireless communication networks: Vision, enabling technologies, and new paradigm shifts,” Sci. China Inf. Sci., vol. 64, no. 1, Jan. 2021.
[3] C.-X. Wang, J. Bian, J. Sun, W. Zhang, and M. Zhang, “A survey of 5G channel measurements and models,” IEEE Commun. Surveys Tuts., vol. 20, no. 4, pp. 3142–3168, 4th Quart., 2018.
[4] P. Qarabaqi and M. Stojanovic, “Statistical characterization and computationally efficient modeling of a class of underwater acoustic communication channels,” IEEE J. Oceanic Eng., vol. 38, no. 4, pp. 701–717, Oct. 2013.
[5] E. Baktash, M. J. Dehghani, M. R. F. Nasab, and M. Karimi, “Shallow water acoustic channel modeling based on analytical second order statistics for moving transmitter/receiver,” IEEE Trans. Signal Process., vol. 63, no. 10, pp. 2533–2545, May 2015.
[6] J. Rudander, T. Husøy, P. Orten, and P. van Walree, “Shallow-water channel sounding for high speed acoustic communication,” in Proc. OCEANS’17, Aberdeen, United Kingdom: IEEE, June 2017, pp. 1–8.
[7] M. Stojanovic and J. Preisig, “Underwater acoustic communication channels: Propagation models and statistical characterization,” IEEE Commun. Mag., vol. 47, no. 1, pp. 84–89, Jan. 2009.
[8] P. A. van Walree, “Propagation and scattering effects in underwater acoustic communication channels,” IEEE J. Oceanic Eng., vol. 38, no. 4, pp. 614–631, Oct. 2013.
[9] P. A. van Walree and R. Otnes, “Ultrawideband underwater acoustic communication channels,” IEEE J. Oceanic Eng., vol. 38, no. 4, pp. 678–688, Oct. 2013.
[10] S. Gul, S. S. H. Zaidi, R. Khan and A. B. Wala, “Underwater acoustic channel modeling using BELLHOP ray tracing method,” in Proc. IBCAST’17, Islamabad, Pakistan, Jan. 2017, pp. 665-670.
[11] J. Huang and R. Diamant, “Pre-setting of channel types for long range underwater acoustic communications,” in Proc. OCEANS’19, Marseille, France, Jun. 2019, pp. 1–6.
[12] F.-X. Socheleau, C. Laot, and J.-M. Passerieux, “A maximum entropy framework for statistical modeling of underwater acoustic communication channels,” in Proc. OCEANS’10, Sydney, Australia: IEEE, May 2010, pp. 1–7.
[13] M. Naderi, M. Pätzold, R. Hicheri, and N. Youssef, “A geometry-based underwater acoustic channel model allowing for sloped ocean bottom conditions,” IEEE Trans. Wireless Commun., vol. 16, no. 4, pp. 2394–2408, Apr. 2017.
[14] M. Naderi, A. G. Zajić, and M. Pätzold, “A nonisovelocity geometry-based underwater acoustic channel model,” IEEE Trans. Veh. Technol., vol. 67, no. 4, pp. 2864–2879, Apr. 2018.
[15] A. G. Zajić, “Statistical modeling of MIMO mobile-to-mobile underwater channels,” IEEE Trans. Veh. Technol., vol. 60, no. 4, pp. 1337–1351, May 2011.
[16] A. Song, M. Stojanovic, and M. Chitre,”Editorial underwater acoustic communications: where we stand and what is next?,” IEEE J. Oceanic Eng., vol. 44, no. 1, pp. 1–6, Jan. 2019.
[17] S. Zhou and Z. Wang, OFDM for Underwater Acoustic Communications. Chichester, West Sussex, United Kingdom: Wiley, 2014.
[18] J. Bian, C.-X. Wang, X. Gao, X. You, and M. Zhang, “A general 3D non-stationary wireless channel model for 5G and beyond,” IEEE Trans. Wireless Commun., accepted for publication.
[19] L. M. Brekhovskikh and Y. P. Lysanov, Fundamentals of Ocean Acoustics, 3rd ed. New York: Springer, 2003.
[20] S. Wu, C.-X. Wang, H. Aggoune, M. M. Alwakeel, and X. You, “A general 3D non-stationary 5G wireless channel model,” IEEE Trans. Commun., vol. 66, no. 7, pp. 3065–3078, July 2018.