RF

    radio frequency
THz

    terahertz
HRTRN

    hybrid RF and THz relay network
5G-NR

    5G New Radio
3GPP

    3rd Generation Partnership Project
ABS

    aerial base station
AC

    address coding
ACF

    autocorrelation function
ACR

    autocorrelation receiver
ADC

    analog-to-digital converter
AIC

    Analog-to-Information Converter
AIC

    Akaike information criterion
ARIC

    asymmetric restricted isometry constant
ARIP

    asymmetric restricted isometry property
ARQ

    Automatic Repeat Request
AUB

    asymptotic union bound
AWGN

    Additive White Gaussian Noise
AWGN

    additive white Gaussian noise
PSK

    asymmetric PSK
AWRICs

    asymmetric weak restricted isometry constants
AWRIP

    asymmetric weak restricted isometry property
BCH

    Bose, Chaudhuri, and Hocquenghem
BCHSC

    BCH based source coding
BEP

    bit error probability
BFC

    block fading channel
BG

    Bernoulli-Gaussian
BGG

    Bernoulli-Generalized Gaussian
BPAM

    binary pulse amplitude modulation
BPDN

    Basis Pursuit Denoising
BPPM

    binary pulse position modulation
BPSK

    Binary Phase Shift Keying
BPZF

    bandpass zonal filter
BSC

    binary symmetric channels
BU

    Bernoulli-uniform
BER

    bit error rate
BS

    base station
BW

    BandWidth
BLLL

     binary log-linear learning
CP

    Cyclic Prefix
CDF

    cumulative distribution function
CDF

    Cumulative Distribution Function
CDF

    cumulative distribution function
CCDF

    complementary cumulative distribution function
CCDF

    complementary CDF
CCDF

    complementary cumulative distribution function
CD

    cooperative diversity
CDMA

    Code Division Multiple Access
ch.f.

    characteristic function
CIR

    channel impulse response
CoSaMP

    compressive sampling matching pursuit
CR

    cognitive radio
CS

    compressed sensing
CS

    Compressed sensing
CS

    compressed sensing
CSI

    channel state information
CCSDS

    consultative committee for space data systems
CC

    convolutional coding
COVID-19

    Coronavirus disease
DAA

    detect and avoid
DAB

    digital audio broadcasting
DCT

    discrete cosine transform
DFT

    discrete Fourier transform
DR

    distortion-rate
DS

    direct sequence
DS-SS

    direct-sequence spread-spectrum
DTR

    differential transmitted-reference
DVB-H

    digital video broadcasting – handheld
DVB-T

    digital video broadcasting – terrestrial
DL

    DownLink
DSSS

    Direct Sequence Spread Spectrum
DFT-s-OFDM

    Discrete Fourier Transform-spread-Orthogonal Frequency Division Multiplexing
DAS

    Distributed Antenna System
DNA

    DeoxyriboNucleic Acid
EC

    European Commission
EED

    exact eigenvalues distribution
EIRP

    Equivalent Isotropically Radiated Power
ELP

    equivalent low-pass
eMBB

    Enhanced Mobile Broadband
EMF

    ElectroMagnetic Field
EU

    European union
EI

    Exposure Index
eICIC

    enhanced Inter-Cell Interference Coordination
FC

    fusion center
FCC

    Federal Communications Commission
FEC

    forward error correction
FFT

    fast Fourier transform
FH

    frequency-hopping
FH-SS

    frequency-hopping spread-spectrum
FS

    Frame synchronization
FS

    frame synchronization
FDMA

    Frequency Division Multiple Access
GA

    Gaussian approximation
GF

    Galois field
GG

    Generalized-Gaussian
GIC

    generalized information criterion
GLRT

    generalized likelihood ratio test
GPS

    Global Positioning System
GMSK

    Gaussian Minimum Shift Keying
GSMA

    Global System for Mobile communications Association
GS

    ground station
GMG

     Grid-connected MicroGeneration
HAP

    high altitude platform
HetNet

    Heterogeneous network
IDR

    information distortion-rate
IFFT

    inverse fast Fourier transform
IHT

    iterative hard thresholding
i.i.d.

    independent, identically distributed
IoT

    Internet of Things
IR

    impulse radio
LRIC

    lower restricted isometry constant
LRICt

    lower restricted isometry constant threshold
ISI

    intersymbol interference
ITU

    International Telecommunication Union
ICNIRP

    International Commission on Non-Ionizing Radiation Protection
IEEE

    Institute of Electrical and Electronics Engineers
ICES

    IEEE international committee on electromagnetic safety
IEC

    International Electrotechnical Commission
IARC

    International Agency on Research on Cancer
IS-95

    Interim Standard 95
KPI

    Key Performance Indicator
LEO

    low earth orbit
LF

    likelihood function
LLF

    log-likelihood function
LLR

    log-likelihood ratio
LLRT

    log-likelihood ratio test
LoS

    Line-of-Sight
LRT

    likelihood ratio test
LWRIC

    lower weak restricted isometry constant
LWRICt

    LWRIC threshold
LPWAN

    Low Power Wide Area Network
LoRaWAN

    Low power long Range Wide Area Network
NLoS

    Non-Line-of-Sight
Li-Fi

    light-fidelity
LED

    light emitting diode
LABS

    LoS transmission with each ABS
NLABS

    NLoS transmission with each ABS
MB

    multiband
MC

    macro cell
MDS

    mixed distributed source
MF

    matched filter
m.g.f.

    moment generating function
MI

    mutual information
MIMO

    Multiple-Input Multiple-Output
MISO

    multiple-input single-output
MJSO

    maximum joint support cardinality
ML

    maximum likelihood
MMSE

    minimum mean-square error
MMV

    multiple measurement vectors
MOS

    model order selection
${M}$-PSK

     $M$-ary phase shift keying
${M}$-PSK

     $M$-ary asymmetric PSK
MP

     multi-period
MINLP

    mixed integer non-linear programming
$M$-QAM

     $M$-ary quadrature amplitude modulation
MRC

    maximal ratio combiner
MSO

    maximum sparsity order
M2M

    Machine-to-Machine
MUI

    multi-user interference
mMTC

    massive Machine Type Communications
mm-Wave

    millimeter-wave
MP

    mobile phone
MPE

    maximum permissible exposure
MAC

    media access control
NB

    narrowband
NBI

    narrowband interference
NLA

    nonlinear sparse approximation
NLOS

    Non-Line of Sight
NTIA

    National Telecommunications and Information Administration
NTP

    National Toxicology Program
NHS

    National Health Service
LOS

    Line of Sight
OC

    optimum combining
OC

    optimum combining
ODE

    operational distortion-energy
ODR

    operational distortion-rate
OFDM

    Orthogonal Frequency-Division Multiplexing
OMP

    orthogonal matching pursuit
OSMP

    orthogonal subspace matching pursuit
OQAM

    offset quadrature amplitude modulation
OQPSK

    offset QPSK
OFDMA

    Orthogonal Frequency-division Multiple Access
OPEX

    Operating Expenditures
OQPSK/PM

    OQPSK with phase modulation
PAM

    pulse amplitude modulation
PAR

    peak-to-average ratio
PDF

    probability density function
PDF

    probability density function
PDF

    probability distribution function
PDP

    power dispersion profile
PMF

    probability mass function
PMF

    probability mass function
PN

    pseudo-noise
PPM

    pulse position modulation
PRake

    Partial Rake
PSD

    power spectral density
PSEP

    pairwise synchronization error probability
PSK

    phase shift keying
PD

    power density
$8$-PSK

     $8$-phase shift keying
PPP

    Poisson point process
PCP

    Poisson cluster process
FSK

    Frequency Shift Keying
QAM

    Quadrature Amplitude Modulation
QPSK

    Quadrature Phase Shift Keying
OQPSK/PM

    OQPSK with phase modulator
RD

    raw data
RDL

    ”random data limit”
RIC

    restricted isometry constant
RICt

    restricted isometry constant threshold
RIP

    restricted isometry property
ROC

    receiver operating characteristic
RQ

    Raleigh quotient
RS

    Reed-Solomon
RSSC

    RS based source coding
r.v.

    random variable
R.V.

    random vector
RMS

    root mean square
RFR

    radiofrequency radiation
RIS

    Reconfigurable Intelligent Surface
RNA

    RiboNucleic Acid
RRM

    Radio Resource Management
RUE

    reference user equipments
RAT

    radio access technology
RB

    resource block
SA-Music

    subspace-augmented MUSIC with OSMP
SC

    small cell
SCBSES

    Source Compression Based Syndrome Encoding Scheme
SCM

    sample covariance matrix
SEP

    symbol error probability
SG

    sparse-land Gaussian model
SIMO

    single-input multiple-output
SINR

    signal-to-interference plus noise ratio
SIR

    signal-to-interference ratio
SISO

    Single-Input Single-Output
SMV

    single measurement vector
SNR

    signal-to-noise ratio
SP

    subspace pursuit
SS

    spread spectrum
SW

    sync word
SAR

    specific absorption rate
SSB

    synchronization signal block
SR

    shrink and realign
tUAV

    tethered Unmanned Aerial Vehicle
TBS

    terrestrial base station
uUAV

    untethered Unmanned Aerial Vehicle
PDF

    probability density functions
PL

    path-loss
TH

    time-hopping
ToA

    time-of-arrival
TR

    transmitted-reference
TW

    Tracy-Widom
TWDT

    TW Distribution Tail
TCM

    trellis coded modulation
TDD

    Time-Division Duplexing
TDMA

    Time Division Multiple Access
Tx

    average transmit
UAV

    Unmanned Aerial Vehicle
URIC

    upper restricted isometry constant
URICt

    upper restricted isometry constant threshold
UWB

    ultrawide band
UWB

    Ultrawide band
URLLC

    Ultra Reliable Low Latency Communications
UWRIC

    upper weak restricted isometry constant
UWRICt

    UWRIC threshold
UE

    User Equipment
UL

    UpLink
WiM

    weigh-in-motion
WLAN

    wireless local area network
WM

    Wishart matrix
WMAN

    wireless metropolitan area network
WPAN

    wireless personal area network
WRIC

    weak restricted isometry constant
WRICt

    weak restricted isometry constant thresholds
WRIP

    weak restricted isometry property
WSN

    wireless sensor network
WSS

    Wide-Sense Stationary
WHO

    World Health Organization
Wi-Fi

    Wireless Fidelity
SpaSoSEnc

    sparse source syndrome encoding
VLC

    Visible Light Communication
VPN

    Virtual Private Network
RF

    Radio Frequency
FSO

    Free Space Optics
IoST

    Internet of Space Things
GSM

    Global System for Mobile Communications
2G

    Second-generation cellular network
3G

    Third-generation cellular network
4G

    Fourth-generation cellular network
5G

    Fifth-generation cellular network
gNB

    next-generation Node-B Base Station
NR

    New Radio
UMTS

    Universal Mobile Telecommunications Service
LTE

    Long Term Evolution
QoS

    Quality of Service

#  Coverage Analysis of Hybrid RF/THz Networks With Best Relay Selection

Zhengying Lou, Baha Eddine Youcef Belmekki, and Mohamed-Slim Alouini, Fellow,
IEEE The authors are with King Abdullah University of Science and Technology
(KAUST), CEMSE division, Thuwal 23955-6900, Saudi Arabia (e-mail:
<EMAIL_ADDRESS><EMAIL_ADDRESS>slim.alouini@kaust.edu.sa).

###### Abstract

Utilizing terahertz (THz) transmission to enhance coverage has proven various
benefits compared to traditional radio frequency (RF) counterparts. This
letter proposes a dual-hop decode-and-forward (DF) routing protocol in a
hybrid RF and THz relay network named hybrid relay selection (HRS). The
coverage probability of the HRS protocol is derived. The HRS protocol
prioritizes THz relays for higher data rates or short source-destination
distances; and RF relays for lower data rates or large source-destination
distances. The proposed HRS protocol offers nearly the same performance as the
optimal selection protocol, which requires complete instantaneous channel
state information (CSI) of all the nodes.

###### Index Terms:

Stochastic geometry, hybrid RF and THz relay network, relay selection,
coverage probability.

##  I Introduction

With the upcoming bandwidth-hungry applications of the sixth generation (6G)
networks such as virtual reality (VR) and holographic communications,
ubiquitous and ultra-high-speed access are required [1]. Utilizing the
terahertz (THz) band will increase channel bandwidth and transmission capacity
as it is considered one of the key enablers of 6G networks [2]. Compared to
traditional radio frequency (RF) networks, dense THz networks have lower power
consumption, smaller time delay, weaker radiation, better beam directivity,
and higher interference immunity [3]. However, THz signals suffer from deep
fading due to large free space path-loss and water molecule absorption, which
significantly limit the effective communication distance [4]. Hence, relay
communications are used to cope with this limitation.

Most of the literature focuses on coverage performance for a given relay link
[5, 6, 7].  The works in [5, 6] studied the outage probability, which is the
complementary cumulative distribution function of the coverage probability, of
dual-hop THz-THz links, and [6] further optimized the transmission power. In
[7], the authors derived the outage probability and the ergodic capacity of
dual-hop THz-RF links, respectively. However, the above distance-fixed three-
point (source-relay-destination) model is simplified for a dense DF THz
network.  The authors in [8] analyzed the downlink coverage probability of the
hybrid RF and THz network, but relay communication is not considered.

Against this background, THz relay selection strategies are lacking in the
literature. Moreover, since THz base stations are expected to coexist with RF
base stations, a relay selection mechanism involving both THz and RF relays is
necessary especially for ultra-dense networks (UDNs) and integrated access and
backhaul (IAB) deployments. In this treatise, we propose a hybrid relay
selection (HRS) protocol with RF and THz relays. Furthermore, the mathematical
derivations are based on a realistic THz channel model that follows
$\alpha-\mu$ which is an experimentally validated model [9], and the relays
are randomly located following a Poisson point process (PPP). The aim of this
letter is to investigate the performance of the proposed protocol and how the
relay selection is carried out when both RF and THz relays are involved in the
selection process. Our proposed HRS protocol is compared with the optimal
relay selection that requires a full and perfect instantaneous channel state
information (CSI) of all nodes in the networks which is impractical to
estimate. We also show when to use RF relays and/or THz relays. Finally, for
the sake of completeness, we compare our protocol with RF and THz direct
transmissions.

##  II System Model

Figure 1: System model.

###  II-A Network Model

We consider the system model depicted in Fig. 1. Without loss of generality,
we set the position of a destination ($D$) at the origin for tractability
since the same performance is obtained for any other $D$ locations due to
Slivnyak’s Theorem. Also, the direction from a source ($S$) to $D$ is the
positive x-axis. In polar coordinate $\left(\rho,\theta\right)$, the positions
of $S$ and $D$ can be denoted as $x_{S}(R_{SD},0)$ and $x_{D}(0,0)$,
respectively, where $R_{SD}$ is the distance between $S$ and $D$.  The
locations of $Q$ relays form a homogeneous PPP denoted by
$\Phi_{Q}=\\{x_{Q,1},x_{Q,2},...\\}$ and with density $\lambda_{Q}$, where $Q$
is replaced by RF in RF networks, and by THz in THz networks.  We assume
$\Phi_{\rm{RF}}$ and $\Phi_{\rm{THz}}$ are independent.  All RF and THz
relays, referred to as all nodes, are located on a two-dimensional circular
disc $\mathcal{C}(x_{D},R_{C})$ centered at the origin with radius $R_{C}$.
Note that the relays outside $\mathcal{C}(x_{D},R_{C})$ are not considered
since they are far from $D$ and suffer a large signal attenuation. Hence,
$\mathcal{C}(x_{D},R_{C})$ is equivalent to the entire ${\mathbb{R}}^{2}$. We
assume that the locations of $S$, $D$, and the relays and the frequency band
of the relays are shared, so the beams have been trained, and the
communication targets are always within the main lobe of beams.  Owing to the
utilization of orthogonal sub-bands and narrow-beam transmission, it is
extremely uncommon that the interference signal in the same sub-band is placed
in the main lobe of the receiving device, hence the interference is minor
compared to noise.

###  II-B Channel Model

####  II-B1 RF Channel Model

RF channels experience path-loss and small-scale fading; therefore, the
corresponding received power can be expressed as,
$h_{\rm{RF}}\left(R_{X}\right)=\varepsilon_{{\rm{RF}}}\,G_{{\rm{RF}}}\,\gamma_{\rm{RF}}\,R_{X}^{-\beta_{\rm{RF}}}\,\mathcal{X}_{\rm{RF}}$,
where $X=\\{S,D\\}$, $R_{X}$ is the distance between a relay and $X$,
$\varepsilon_{{\rm{RF}}}$ is the transmitting power, and $G_{{\rm{RF}}}$ is
total antenna gain. The path-loss is modeled as
$\gamma_{\rm{RF}}R_{X}^{-\beta_{\rm{RF}}}$, where $\beta_{\rm{RF}}$ is the
path-loss exponent,
$\gamma_{\rm{RF}}=\left({c}/{4\pi\nu_{\rm{RF}}}\right)^{2}$, $c=3\times
10^{8}$ m/s is the speed of light, and $\nu_{\rm{RF}}$ is the RF carrier
frequency. The small-scale fading $\mathcal{X}_{\rm{RF}}$ is subject to an
exponential distribution with unit mean. The signal-to-noise ratio (SNR) of
the RF channel is given by

$\small{\rm{SNR}}_{{\rm{RF}},X}=\frac{\varepsilon_{{\rm{RF}}}\,G_{\rm{RF}}\gamma_{\rm{RF}}\,R_{X}^{-\beta_{\rm{RF}}}\,\mathcal{X}_{\rm{RF}}}{\sigma_{\rm{RF}}^{2}},\vspace{-1mm}$
(1)

where $\sigma_{\rm{RF}}^{2}$ is the thermal noise. We note that
$\sigma_{\rm{RF}}^{2}$ is a function of the transmission bandwidth
$B_{\rm{RF}}$.

####  II-B2 THz Channel Model

According to [8], the received power in THz transmission is modeled as,
$h_{\rm{THz}}\left(R_{X}\right)=\varepsilon_{{\rm{THz}}}\,G_{{\rm{THz}}}\,\gamma_{{\rm{THz}}}\,\mathcal{X}_{\rm{THz}}\,{\exp\left({-\beta_{\rm{THz}}R_{X}}\right)}/{R_{X}^{2}}$,
where $X=\\{S,D\\}$, $R_{X}$ is the distance between a relay and $X$,
$\varepsilon_{{\rm{THz}}}$ is the transmitting power, and $G_{{\rm{THz}}}$ is
the total antenna gain. The path-loss is modeled as
$\gamma_{{\rm{THz}}}{\exp\left({-\beta_{\rm{THz}}R_{X}}\right)}/{R_{X}^{2}}$,
where $\gamma_{{\rm{THz}}}=\left({c}/{4\pi\nu_{\rm{THz}}}\right)^{2}$ ,
$\nu_{\rm{THz}}$ is the THz carrier frequency in GHz, and $\beta_{\rm{THz}}$
is the molecular absorption coefficient related to $\nu_{\rm{THz}}$. The
small-scale fading $\mathcal{X}_{\rm{THz}}$ is modeled as $\alpha-\mu$
distribution which is an experimentally validated fading model for THz
frequencies [9]. Moreover, $\alpha$ denotes the fading parameter and $\mu$
denotes the normalized variance of the channel fading. The complementary
cumulative distribution function (CCDF) of $\mathcal{X}_{\rm{THz}}$ is given
by $\bar{F}_{\mathcal{X}_{\rm{THz}}}(m)={\Gamma(\mu,\mu
m^{\frac{\alpha}{2}})}/{\Gamma(\mu)}$, where ${\Gamma(\mu,\mu
m^{\frac{\alpha}{2}})}$ and ${\Gamma(\mu)}$ are the upper incomplete Gamma
function and the Gamma function, respectively [10]. The SNR of the THz channel
is given by

${\rm{SNR}}_{{\rm{THz}},X}=\frac{\varepsilon_{{\rm{THz}}}\,G_{{\rm{THz}}}\gamma_{{\rm{THz}}}\,\exp\left({-\beta_{\rm{THz}}R_{X}}\right)\,\mathcal{X}_{\rm{THz}}}{\sigma_{\rm{THz}}^{2}\,R_{X}^{2}},\vspace{-1mm}$
(2)

where $\sigma_{\rm{THz}}^{2}$ is the thermal noise, which is a function of the
transmission bandwidth $B_{\rm{THz}}$.

###  II-C Relay Selection Strategy

Due to excessive path-loss and blockage, a half-duplex DF relay protocol is
adopted since the direct link may be blocked. In fixed topology DF networks,
the optimal relay selection relies on the maximization of end-to-end
${\rm{SNR}}_{Q,X}$ or the maximization of the minimum ${\rm{SNR}}_{Q,X}$ of
$S-x_{Q,i}$ and $x_{Q,i}-D$ links. However, this protocol requires a full and
perfect instantaneous CSI of all the nodes in the networks, that is, CSI
between $S$ and all relays, and between all relays and $D$. In addition, due
to the randomness of the relay locations, analysis is intractable [11].
Therefore, we proposed a slightly suboptimal protocol that only requires
instantaneous CSI between $S$ and all nodes while still being mathematically
tractable.  The protocol has the following selection strategy: we first select
two sets of relays $\widetilde{\Phi}_{Q}$ that can meet the attempted rate
$y_{th}$ of $S$. Among these sets, only the relay
$x_{\rm{RF},i}\in\widetilde{\Phi}_{\rm{RF}}$ or
$x_{\rm{THz},i}\in\widetilde{\Phi}_{\rm{THz}}$ that has the best link to $D$,
that is, provides the maximum average rate with $D$ is selected [8]. The set
of relays that can establish reliable communication with $S$ can be expressed
as

$\widetilde{\Phi}_{Q}=\\{x_{Q,i}\in\Phi_{Q},\,{\rm{SNR}}_{Q,S}\geq\tau_{Q}\\},\
Q=\\{\rm{RF},\rm{THz}\\},\vspace{-1mm}$ (3)

where $\tau_{Q}$ is a predefined threshold, ${\rm{SNR}}_{Q,S}$ is the SNR of
S-relay link. Subsequently, the relation between achievable rate $y_{Q,X}$ and
${\rm{SNR}}_{Q,X}$ is given by

$y_{Q,X}=\frac{B_{Q}}{2}\log_{2}\left(1+{\rm{SNR}}_{Q,X}\right),\vspace{-1mm}$
(4)

where $B_{Q}$ is the transmission bandwidth of $Q$ link and the factor 1/2 is
used because two time slots are required to transmit data via the relay node.
The rate coverage probability is defined as the probability of achieving the
desired rate $y_{th}$, i.e.,
${\textrm{SNR}}_{Q,X}>\tau_{Q}=2^{\left(2y_{th}/B_{Q}\right)}-1$.

##  III Coverage Probability of HRS Protocol

In this section, we derive the coverage probability, the intermediate distance
distributions, and the association probabilities of the HRS protocol. Note
that the nearest relay provides the largest average rate for $D$ in our model,
so the distance distributions are given in Lemmas.

######  Lemma 1.

The complementary cumulative distribution function (CCDF) of the distance
between D and its nearest RF relay is given by

$\small\begin{split}\overline{F}_{R_{D},{\rm{RF}}}(r)=\exp\bigg{(}-\int_{0}^{r}\int_{0}^{2\pi}\lambda_{\rm{RF}}\exp\Big{(}-\tau_{\rm{RF}}\varepsilon_{{\rm{RF}}}^{-1}\\\
\times
G_{{\rm{RF}}}^{-1}\gamma_{{\rm{RF}}}^{-1}\\!\left(\rho^{2}\\!+\\!R_{SD}^{2}\\!-\\!2\rho
R_{SD}\cos\theta\right)^{\frac{\beta_{\rm{RF}}}{2}}\sigma_{\rm{RF}}^{2}\Big{)}\bigg{)}\rho\,\mathrm{d}\theta\mathrm{d}\rho.\end{split}\vspace{-1mm}$
(5)

where $0<r\leq R_{C}$.

Then, the probability density functions (PDF) of the distance between D and
its nearest RF relay is given by,

$\small\begin{split}&f_{R_{D},{\rm{RF}}}(r)=r\lambda_{\rm{RF}}\int_{0}^{2\pi}\exp\Big{(}-\tau_{\rm{RF}}\varepsilon_{{\rm{RF}}}^{-1}G_{{\rm{RF}}}^{-1}\gamma_{{\rm{RF}}}^{-1}\\\
&\times\left(\rho^{2}+R_{SD}^{2}-2\rho
R_{SD}\cos\theta\right)^{\frac{\beta_{\rm{RF}}}{2}}\sigma_{\rm{RF}}^{2}\Big{)}\mathrm{d}\theta\\\
&\times\exp\bigg{(}-\int_{0}^{r}\int_{0}^{2\pi}\lambda_{\rm{RF}}\exp\Big{(}-\tau_{\rm{RF}}\varepsilon_{{\rm{RF}}}^{-1}G_{{\rm{RF}}}^{-1}\gamma_{{\rm{RF}}}^{-1}\\\
&\times\left(\rho^{2}+R_{SD}^{2}-2\rho
R_{SD}\cos\theta\right)^{\frac{\beta_{\rm{RF}}}{2}}\sigma_{\rm{RF}}^{2}\Big{)}\bigg{)}\rho\,\mathrm{d}\theta\mathrm{d}\rho.\\\
\end{split}\vspace{-1mm}$ (6)

###### Proof.

The CCDF of the distance can be derived by using the void probability of a
homogeneous PPP [12]

$\small\begin{split}\overline{F}_{R_{D},{\rm{RF}}}(r)=\mathbb{P}\left[{\rm
No\,relays\,in\,}\mathcal{C}(x_{D},r)\right]=\exp\left(-\Lambda_{\rm{RF}}\left(r\right)\right),\end{split}$
(7)

where $\mathcal{C}(x_{D},r)$ is the disk centered at the origin ($x_{D}$) with
radius $r$. The mean number of RF relays
$x_{\rm{RF},i}\in\widetilde{\Phi}_{\rm{RF}}$ in the circular disc
$\mathcal{C}(x_{D},r)$, denoted by $\Lambda_{\rm{RF}}\left(r\right)$, can be
written as

$\Lambda_{\rm{RF}}\left(r\right)=\int_{0}^{r}\int_{0}^{2\pi}\widehat{\lambda}_{\rm{RF}}(\rho,\theta)\rho\,\mathrm{d}\theta\mathrm{d}\rho,\vspace{-1mm}$
(8)

where $\rho$ and $\theta$ are the radial distance and polar angle in the polar
coordinate system respectively, and $\widehat{\lambda}_{\rm{RF}}(\rho,\theta)$
is the density of $\widetilde{\Phi}_{\rm{RF}}$. The set
$\widetilde{\Phi}_{\rm{RF}}$ can be regarded as a dependent-thinning of
$\Phi_{\rm{RF}}$, and only the relays that have an ${\rm{SNR}}_{{\rm{RF}},S}$
greater than the threshold $\tau_{\rm{RF}}$ are retained. Hence, the set
$\widetilde{\Phi}_{\rm{RF}}$ is a inhomogeneous PPP of density
$\widehat{\lambda}_{\rm{RF}}(\rho,\theta)$. The density
$\widehat{\lambda}_{\rm{RF}}(\rho,\theta)$ can be derived as follows

$\small\begin{split}\widehat{\lambda}_{\rm{RF}}(\rho,\theta)&=\lambda_{\rm{RF}}\,\mathbb{P}\left[\frac{\varepsilon_{{\rm{RF}}}\,G_{\rm{RF}}\gamma_{\rm{RF}}\,R_{S}^{-\beta_{\rm{RF}}}\,\mathcal{X}_{\rm{RF}}}{\sigma_{\rm{RF}}^{2}}>\tau_{\rm{RF}}\right]\\\
&=\lambda_{\rm{RF}}\,\mathbb{P}\left[\mathcal{X}_{\rm{RF}}>\tau_{\rm{RF}}\varepsilon_{{\rm{RF}}}^{-1}G_{{\rm{RF}}}^{-1}\gamma_{{\rm{RF}}}^{-1}R_{S}^{\beta_{\rm{RF}}}\sigma_{\rm{RF}}^{2}\right]\\\
&=\lambda_{\rm{RF}}\exp\left(-\tau_{\rm{RF}}\varepsilon_{{\rm{RF}}}^{-1}G_{{\rm{RF}}}^{-1}\gamma_{{\rm{RF}}}^{-1}R_{S}^{\beta_{\rm{RF}}}\sigma_{\rm{RF}}^{2}\right),\end{split}\vspace{-2mm}$
(9)

where $R_{S}$ is the distance between $S$ and relay and can expressed as

$R_{S}^{2}=\rho^{2}+R_{SD}^{2}-2\rho\,R_{SD}\cos\theta.\vspace{-1mm}$ (10)

By substituting (8), (9), and (10) into (7), we obtain (5). Then, using
$f_{R_{D},{\rm{RF}}}(r)=-\frac{\mathrm{d}}{\mathrm{d}r}\overline{F}_{R_{D},{\rm{RF}}}(r)$
and Leibniz integral rule, (6) is derived straightforwardly.

∎

######  Lemma 2.

The CCDF and PDF of the distance between D and its nearest THz relay are given
by (11) and (12), respectively, shown at the top of the next page.

$\small\begin{split}\overline{F}_{R_{D},{\rm{THz}}}(r)=\exp\left(-\frac{\lambda_{\rm{THz}}}{\Gamma\left(\mu\right)}\int_{0}^{r}\int_{0}^{2\pi}{\Gamma\left(\mu,\mu\,\left(\frac{\tau_{\rm{THz}}\left(\rho^{2}+R_{SD}^{2}-2\rho
R_{SD}\cos\theta\right)\sigma_{\rm{THz}}^{2}}{\varepsilon_{{\rm{THz}}}G_{{\rm{THz}}}\gamma_{{\rm{THz}}}\exp\left(-\beta_{\rm{THz}}\left(\rho^{2}+R_{SD}^{2}-2\rho
R_{SD}\cos\theta\right)^{\frac{1}{2}}\right)}\right)^{\frac{\alpha}{2}}\right)}\rho\,\mathrm{d}\theta\mathrm{d}\rho\right).\end{split}$
(11)

$\small\begin{split}f_{R_{D},{\rm{THz}}}(r)&=r\,\frac{\lambda_{\rm{THz}}}{\Gamma\left(\mu\right)}\int_{0}^{2\pi}{\Gamma\left(\mu,\mu\,\left(\frac{\tau_{\rm{THz}}\left(r^{2}+R_{SD}^{2}-2rR_{SD}\cos\theta\right)\sigma_{\rm{THz}}^{2}}{\varepsilon_{{\rm{THz}}}G_{{\rm{THz}}}\gamma_{{\rm{THz}}}\exp\left(-\beta_{\rm{THz}}\left(r^{2}+R_{SD}^{2}-2rR_{SD}\cos\theta\right)^{\frac{1}{2}}\right)}\right)^{\frac{\alpha}{2}}\right)}\,\mathrm{d}\theta\\\
\vspace{-1mm}&\times\exp\left(-\frac{\lambda_{\rm{THz}}}{\Gamma\left(\mu\right)}\int_{0}^{r}\int_{0}^{2\pi}{\Gamma\left(\mu,\mu\,\left(\frac{\tau_{\rm{THz}}\left(\rho^{2}+R_{SD}^{2}-2\rho
R_{SD}\cos\theta\right)\sigma_{\rm{THz}}^{2}}{\varepsilon_{{\rm{THz}}}G_{{\rm{THz}}}\gamma_{{\rm{THz}}}\exp\left(-\beta_{\rm{THz}}\left(\rho^{2}+R_{SD}^{2}-2\rho
R_{SD}\cos\theta\right)^{\frac{1}{2}}\right)}\right)^{\frac{\alpha}{2}}\right)}\rho\,\mathrm{d}\theta\mathrm{d}\rho\right).\end{split}\vspace{-0.2cm}$
(12)

###### Proof.

Following the same steps as in Lemma 1, we obtain

$\small\overline{F}_{R_{D},{\rm{THz}}}(r)=\exp\left({-\int_{0}^{r}\int_{0}^{2\pi}\widehat{\lambda}_{\rm{THz}}(\rho,\theta)\rho\,\mathrm{d}\theta\mathrm{d}\rho}\right),\vspace{-1mm}$
(13)

where $\widehat{\lambda}_{\rm{THz}}(\rho,\theta)$ is the density of
$\widetilde{\Phi}_{\rm{THz}}$, which is given by

$\small\begin{split}&\widehat{\lambda}_{\rm{THz}}(\rho,\theta)=\lambda_{\rm{THz}}\,\mathbb{P}\left[\frac{\varepsilon_{{\rm{THz}}}\,G_{{\rm{THz}}}\gamma_{{\rm{THz}}}\,\,\mathcal{X}_{\rm{THz}}}{\sigma_{\rm{THz}}^{2}\,\exp\left({\beta_{\rm{THz}}R_{S}}\right)\,R_{S}^{2}}>\tau_{\rm{THz}}\right]\\\
&=\lambda_{\rm{THz}}\,\mathbb{P}\bigg{[}\mathcal{X}_{\rm{THz}}>\frac{\tau_{\rm{THz}}R_{S}^{2}\exp\left(\beta_{\rm{THz}}R_{S}\right)\sigma_{\rm{THz}}^{2}}{\varepsilon_{{\rm{THz}}}G_{{\rm{THz}}}\gamma_{{\rm{THz}}}}\bigg{]}\\\
&\overset{(a)}{=}\frac{\lambda_{\rm{THz}}}{\Gamma\left(\mu\right)}\Gamma\left(\mu,\mu\,\left(\frac{\tau_{\rm{THz}}R_{S}^{2}\exp\left(\beta_{\rm{THz}}R_{S}\right)\sigma_{\rm{THz}}^{2}}{\varepsilon_{{\rm{THz}}}G_{{\rm{THz}}}\gamma_{{\rm{THz}}}}\right)^{\frac{\alpha}{2}}\right),\end{split}\vspace{-1mm}$
(14)

where (a) follows the definition of the CCDF of $\mathcal{X}_{\rm{THz}}$.
Substituting (14) and (10) into (13) concludes the proof of (11). Similar to
the proof of (6) in Lemma 1, (11) is given by
$f_{R_{D},{\rm{THz}}}(r)=-\frac{\mathrm{d}}{\mathrm{d}r}\overline{F}_{R_{D},{\rm{THz}}}(r)$.
∎

Using the above Lemmas, we derive the distribution of RF and THz relays that
are closest to $D$. However, for the nearest RF relay to be selected, it is
necessary to ensure that the achievable rate of THz transmission is less than
that of RF transmission, and vice versa. Therefore, we calculate the
(selection) association probabilities in the following two Lemmas.

######  Lemma 3.

Given that the distance between $D$ and its nearest RF relay is $r$, the
probability that $D$ selects/associates with this RF relay is given by

$P_{\rm{RF}}^{A}(r)=1-\int_{0}^{R_{\rm{R2T}}(r)}f_{R_{D},{\rm{THz}}}(\rho)\mathrm{d}\rho,\vspace{-1mm}$
(15)

where $R_{\rm{R2T}}$ is given by

$\small\begin{split}&R_{\rm{R2T}}(r)=\frac{2}{\beta_{\rm{THz}}}\times\\\
&\mathcal{W}\\!\\!\left(\\!\\!\frac{\beta_{\rm{THz}}}{2}\\!\\!\left(\\!\\!\\!\left(\\!\\!\left(\\!1\\!\\!+\\!\frac{\varepsilon_{{\rm{RF}}}G_{{\rm{RF}}}\gamma_{{\rm{RF}}}}{r^{\beta_{\rm{RF}}}\sigma_{{\rm{RF}}}^{2}}\\!\right)^{\\!\\!\frac{B_{{\rm{RF}}}}{B_{{\rm{THz}}}}}\\!\\!\\!\\!\\!-\\!1\\!\\!\right)\\!\\!\frac{\sigma_{{\rm{THz}}}^{2}}{\varepsilon_{{\rm{THz}}}G_{{\rm{THz}}}\gamma_{{\rm{THz}}}}\\!\\!\right)^{\\!\\!\\!-\\!\frac{1}{2}}\\!\right)\\!.\end{split}\vspace{-1mm}$
(16)

We note that $\mathcal{W}\left(\cdot\right)$ is the Lambert
$\mathcal{W}$-function defined as the inverse of the function
$f_{\mathcal{W}}(m)=me^{m}$.

###### Proof.

The distance of the nearest THz relay, denoted by $R_{\rm{R2T}}$, can be
solved by setting the average achievable rate of the RF relay with distance
$r$ to $D$ equals to that of the THz relay with distance $R_{\rm{R2T}}$ to $D$
[13]

$\begin{split}\mathbb{E}_{\mathcal{X}_{\rm{RF}}}&\left[{\rm{SNR}}_{{\rm{RF}},D}|R_{{\rm{RF}},D}=r\right]\\\
&=\mathbb{E}_{\mathcal{X}_{\rm{THz}}}\left[{\rm{SNR}}_{{\rm{THz}},D}|R_{{\rm{THz}},D}=R_{\rm{R2T}}\right].\end{split}\vspace{-1mm}$
(17)

From the above equation, the relationship between $R_{\rm{R2T}}$ and $r$ is
given by (16). By the definition of association probability
$P_{\rm{RF}}^{A}(r)$ in Lemma 3 and the maximum average received power
associated strategy, the following result is given

$\begin{split}P_{\rm{RF}}^{A}(r)&=\mathbb{P}\left[{\rm
No\,relays\,in\,}\mathcal{C}(x_{D},R_{\rm{R2T}})\right]=\\\
&=1-\int_{0}^{R_{\rm{R2T}}(r)}f_{R_{D},{\rm{THz}}}(\rho)\mathrm{d}\rho.\end{split}\vspace{-1mm}$
(18)

This concludes the proof of Lemma 3. ∎

######  Lemma 4.

Given that the distance between $D$ and its nearest THz relay is $r$, the
probability that $D$ selects/associates with this THz relay is given by

$P_{\rm{THz}}^{A}(r)=1-\int_{0}^{R_{\rm{T2R}}(r)}f_{R_{D},{\rm{RF}}}(\rho)\mathrm{d}\rho,\vspace{-1mm}$
(19)

where $f_{R_{D},{\rm{RF}}}(r)$ is defined in (6) and $R_{\rm{T2R}}(r)$ is
given by

$\small\begin{split}&R_{\rm{T2R}}(r)=\\\
&\left(\\!\\!\\!\left(\\!\\!\\!\left(\\!\\!1\\!\\!+\\!\\!\frac{\varepsilon_{{\rm{THz}}}G_{\rm{THz}}\gamma_{\rm{THz}}}{r^{2}\exp\left(\beta_{\rm{THz}}r\right)\sigma_{\rm{THz}}^{2}}\\!\right)^{\frac{{B_{\rm{THz}}}}{{B_{\rm{RF}}}}}\\!\\!\\!\\!\\!-\\!1\\!\\!\right)\\!\frac{\sigma_{\rm{RF}}^{2}}{\varepsilon_{{\rm{RF}}}G_{\rm{RF}}\gamma_{\rm{RF}}}\\!\\!\right)^{-\frac{1}{\beta_{\rm{RF}}}}\\!.\end{split}\vspace{-1mm}$
(20)

###### Proof.

The proof of Lemma 4 is similar to that of Lemma 3 and therefore, omitted
here. ∎

Figure 2: Density of RF and THz relays under the same coverage probability of
HRS.

Figure 3: Impact of transmission power and desired rate threshold on HRS.

Figure 4: Impact of relay density and source-destination distance on HRS.

Based on the distance distributions and association probabilities, the main
result of this letter is given by the following theorem.

######  Theorem 1.

The coverage probability at $D$ using the HRS protocol is given by

$\small\begin{split}P^{C}=&\int_{0}^{R_{C}}f_{R_{D},{\rm{RF}}}(\rho)\exp\left(-\frac{\tau_{\rm{RF}}\rho^{\beta_{\rm{RF}}}\sigma_{\rm{RF}}^{2}}{\varepsilon_{{\rm{RF}}}G_{{\rm{RF}}}\gamma_{{\rm{RF}}}}\right)P_{\rm{RF}}^{A}\left(\rho\right)\,\mathrm{d}\rho\\\
&+\int_{0}^{R_{C}}f_{R_{D},{\rm{THz}}}(\rho)P_{\rm{THz}}^{A}\left(\rho\right)\frac{1}{\Gamma\left(\mu\right)}\\\
&\times{\Gamma\\!\left(\\!\mu,\mu\,\left(\frac{\tau_{\rm{THz}}\rho^{2}\sigma_{\rm{THz}}^{2}}{\varepsilon_{{\rm{THz}}}G_{{\rm{THz}}}\gamma_{{\rm{THz}}}\exp\left(-\beta_{\rm{THz}}\rho\right)}\right)^{\frac{\alpha}{2}}\right)}\,\mathrm{d}\rho.\end{split}\vspace{-1mm}$
(21)

###### Proof.

The coverage probability $P^{C}$ is defined as the probability that the
received ${\rm{SNR}}_{Q,D}$ at the receiver $D$ from the selected $Q$ relay is
above the predefined threshold $\tau_{Q}$, i.e.,
$P^{C}=\mathbb{P}\left[{\rm{SNR}}_{Q,D}>\tau_{Q}\right]$. We divide the
overall coverage probability $P^{C}$ into coverage probability associated with
RF relay $P_{\rm{RF}}^{C}$ and coverage probability associated with THz relay
$P_{\rm{THz}}^{C}$, i.e., $P^{C}=P_{\rm{RF}}^{C}+P_{\rm{THz}}^{C}$.
Furthermore, $P_{\rm{RF}}^{C}$ can be obtained by (i) solving the coverage
probability given that the distance of $S-x_{\rm{RF},i}$ is
$R_{{\rm{RF}},D}=\rho$; (ii) multiplying the coverage probability with the
association probability given that the distance of $S-x_{\rm{RF},i}$ is
$R_{{\rm{RF}},D}=\rho$ (which is denoted as
$P_{\rm{RF}}^{A}\left(\rho\right)$); (iii) compute the expectation
$\mathbb{E}_{\rho}$. Hence, $P_{\rm{RF}}^{C}$ is calculated as follows

$\small\begin{split}P_{\rm{RF}}^{C}&=\mathbb{E}_{\rho}\left[\mathbb{P}\left[{\rm{SNR}}_{{\rm{RF}},D}>\tau_{\rm{RF}}|R_{{\rm{RF}},D}=\rho\right]P_{\rm{RF}}^{A}\left(\rho\right)\right]\\\
&=\mathbb{E}_{\rho}\left[\exp\left(-\frac{\tau_{\rm{RF}}\rho^{\beta_{\rm{RF}}}\sigma_{\rm{RF}}^{2}}{\varepsilon_{{\rm{RF}}}G_{{\rm{RF}}}\gamma_{{\rm{RF}}}}\right)P_{\rm{RF}}^{A}\left(\rho\right)\right]\\\
&=\int_{0}^{R_{C}}\\!\\!\\!f_{R_{D},{\rm{RF}}}(\rho)\exp\left(-\frac{\tau_{\rm{RF}}\rho^{\beta_{\rm{RF}}}\sigma_{\rm{RF}}^{2}}{\varepsilon_{{\rm{RF}}}G_{{\rm{RF}}}\gamma_{{\rm{RF}}}}\right)P_{\rm{RF}}^{A}\left(\rho\right)\,\mathrm{d}\rho.\end{split}$
(22)

The proof of $P_{\rm{THz}}^{C}$ is similar to the proof of $P_{\rm{RF}}^{C}$,
therefore, it is omitted here.  We note that the association probability
$P_{Q}^{A}$ is contained in the coverage probability associated with $Q$ relay
$P_{Q}^{C}$. By using $P^{C}=P_{\rm{RF}}^{C}+P_{\rm{THz}}^{C}$, the overall
coverage probability is derived. ∎

##  IV Numerical Results

In this section, we provide selected simulation results to demonstrate the
performance of the proposed HRS protocol. The simulations coincide perfectly
with the theoretical analysis, and each simulation is performed over $10^{6}$
independent network realizations within a circular disc of radius $R_{C}=200$
m. We consider the THz link with an antenna gain $G_{\rm{THz}}=40$ dBi and a
carrier frequency $\nu_{\rm{THz}}=1.8$ THz. The corresponding absorption value
$\beta_{\rm{THz}}$ is chosen from the realistic database, which is $0.2$ m-1
[14]. To simulate the $\alpha-\mu$ channel, we set $\alpha=2$ and $\mu=4$ for
the THz link. The carrier frequency $\nu_{\rm{RF}}$ of the RF link is $2.1$
GHz. The antenna gain $G_{\rm{RF}}$ is $20$ dBi and path-loss exponent
$\beta_{\rm{RF}}=2.5$.  We set the noise power to $-174$ dBm/Hz for both links
such as $\sigma^{2}_{Q}=-174+10\log_{10}B_{Q}$ dBm , while the transmission
bandwidths $B_{\rm{RF}}=40$ MHz and $B_{\rm{THz}}=0.5$ GHz, by calculation
$\sigma^{2}_{\rm{RF}}=-98$ dBm and $\sigma^{2}_{\rm{THz}}=-87$ dBm.

Fig. 4 shows the density of RF and THz relays in the HRS protocol required to
achieve a given coverage probability when $\varepsilon_{{Q}}=$ $1$ W,
$R_{SD}=50$ m, and $y_{th}=420$ Mbps. For a given coverage probability, when
$\lambda_{\rm{THz}}$ increases linearly, $\lambda_{\rm{RF}}$ decreases roughly
linearly. Moreover, when the coverage probability grows from $0.7$ to $0.9$,
the slopes of the three curves increase consecutively, but the amplitude of
the increase is quite small. Specifically, approximately $20$ additional THz
relays are required for $1$ RF relay reduced to maintain the same coverage
probability for an area of $2\times 10^{4}$ m2.

(a) $R_{SD}=20$ m (b) $R_{SD}=50$ m (c) $R_{SD}=80$ m

Figure 5: Impact of different transmissions, the source-destination distance
($R_{SD}$), and desired data rate on coverage probability.

In Fig. 4, we keep $R_{SD}=50$ m and select a set of densities with coverage
of 90% as shown in Fig. 4, i.e., $\lambda_{\rm{RF}}=5\times 10^{-4}$ m-2 and
$\lambda_{\rm{THz}}=4\times 10^{-3}$ m-2 in HRS protocol. As shown in Fig. 4,
the coverage probability drops dramatically as the desired rate threshold
$y_{th}$ increases from $400$ Mbps to $800$ Mbps. Moreover, comparing the blue
and purple curves reveals that raising the transmission power of RF relays
only improves performance when the data rate is below $500$ Mbps, but has
essentially no effect when it surpasses $500$ Mbps. In general, increasing the
transmission power of THz relays has a greater impact on performance than
raising that of RF relays, particularly when the data rate is high.

In Fig. 4, we set $\varepsilon_{{Q}}=$ 1 W and select three sets of densities
with 70% coverage as shown in Fig. 4, i.e.,
$\lambda_{\rm{RF}}~{}=~{}\\{0.45,\>1.28,\>2.27\\}\times 10^{-4}$ m-2, and
$\lambda_{\rm{THz}}~{}=~{}\\{5,\;3,\;1\\}\times 10^{-3}$ m-2, respectively. As
is shown in Fig. 4, the HRS coverage probability of THz relay dominated-
network ($\lambda_{\rm{THz}}=5\times 10^{-3}$ m-2) is strongly influenced by
the distance between source and destination $R_{SD}$, especially when
$40\rm{m}\leq R_{SD}\leq 60\rm{m}$. The coverage probability of THz relay-
dominated network decreases drastically as $R_{SD}$ increases compared to RF
relay-dominated network.

In Fig. 5, we compare the proposed HRS protocol with the optimal relay
selection strategy according to the max-min selection scheme. We also compare
the proposed protocol with RF-only relay selection and THz-only relay
selection.  For the sake of fairness, we compare the HRS protocol with RF and
THz direct transmissions, which use the same channel model as the relay link.
Notably, direct transmission takes only one time slot, whereas relay
transmission requires two [15].  We set $\varepsilon_{{Q}}=1$ W and the
density of RF and THz relays are $\lambda_{\rm{RF}}=5\times 10^{-4}$ m-2 and
$\lambda_{\rm{THz}}=4\times 10^{-3}$ m-2. We can notice in Fig. 5(a) that for
short distances ($R_{SD}=20$ m), RF relay selection and RF direct transmission
have lower coverage probability and it starts decreasing drastically for
higher data rates ($500$ Mbps), while $P^{C}\approx 1$ for other four
transmission schemes. Hence, we conclude from Fig. 5(a) that for short
distances: (1) it is better to use THz direct transmission when available than
relay transmission; (2) the HRS protocol offers the same performance as the
optimal protocol without requiring full and perfect CSI. In Fig.5(b), we can
notice that for $R_{SD}=50$ m, the THz direct transmission has zero coverage
probability.  We can also notice that for lower data rates ($300-450$ Mbps),
RF relay and RF direct transmission offer better performance than THz relay.
In contrast, for higher data rates ($500$ Mbps), the HRS protocol will select
only THz relays since THz relay offers better performance since RF relay and
RF direct transmission have a zero coverage probability. Finally, we notice
that the proposed HRS protocol and the optimal protocol have nearly the same
performance; and they both offer better performance than all the
aforementioned schemes. In Fig.5(c), we can see that for large distances
$R_{SD}=80$ m, THz communications (direct and relay) have $P^{C}=0$, and the
HRS protocol will select only RF relays.

##  V Conclusion

In this letter, the coverage probability expression of a proposed hybrid RF
and THz relay selection protocol had been derived. The HRS protocol had mainly
selected THz relays (associated with THz relays) for higher data rates or
short source-destination distances, while it had mainly selected RF relays for
lower data rates or large source-destination distances. The superiority of the
HRS protocol had been demonstrated by comparing it with different strategies.
In addition, the proposed HRS protocol had been compared with the optimal
selection protocol, which requires a full and perfect instantaneous CSI of all
the nodes in the networks. The HRS protocol had offered nearly the same
performance as the optimal selection protocol for lower data rates, and the
same performance for higher data rates.

## References

  * [1] R. Wang, M. A. Kishk, and M.-S. Alouini, “Ultra-dense LEO satellite-based communication systems: A novel modeling technique,” _IEEE Commun. Mag._ , vol. 60, no. 4, pp. 25–31, Apr. 2022.
  * [2] H. Sarieddeen _et al._ , “Next generation Terahertz communications: A rendezvous of sensing, imaging, and localization,” _IEEE Commun. Mag._ , vol. 58, no. 5, pp. 69–75, May 2020.
  * [3] T. S. Rappaport _et al._ , “Millimeter wave mobile communications for 5G cellular: It will work!” _IEEE Access_ , vol. 1, pp. 335–349, May 2013.
  * [4] G. Chattopadhyay _et al._ , “Compact Terahertz instruments for planetary missions,” in _Proc. IEEE 9th Euro. Conf. on Antennas and Propagation_ , Apr. 2015.
  * [5] A.-A. A. Boulogeorgos and A. Alexiou, “Outage probability analysis of THz relaying systems,” in _IEEE 31st Annu. Int. Symp. Personal, Indoor and Mobile Radio Commun._ , Oct. 2020, pp. 1–7.
  * [6] S. Farrag _et al._ , “Outage probability analysis of UAV assisted mobile communications in THz channel,” in _Proc. IEEE 16th Annu. Conf. Wireless On-Demand Netw. Syst. Services Conf._ , Mar. 2021, pp. 1–8.
  * [7] S. Li and L. Yang, “Performance analysis of dual-hop THz transmission systems over $\alpha-\mu$ fading channels with pointing errors,” _IEEE Internet Things J._ , pp. 1–1, Dec. 2021.
  * [8] J. Sayehvand and H. Tabassum, “Interference and coverage analysis in coexisting RF and dense Terahertz wireless networks,” _IEEE Wireless Commun. Lett._ , vol. 9, no. 10, pp. 1738–1742, Oct. 2020.
  * [9] E. N. Papasotiriou _et al._ , “An experimentally validated fading model for thz wireless systems,” _Sci. Rep._ , vol. 11, no. 1, pp. 1–14, Sep. 2021\.
  * [10] H. Lei _et al._ , “Secrecy capacity analysis over $\alpha-\mu$ fading channels,” _IEEE Commun. Lett._ , vol. 21, no. 6, pp. 1445–1448, Jun. 2017\.
  * [11] R. Wang, M. A. Kishk, and M.-S. Alouini, “Stochastic geometry-based low latency routing in massive LEO satellite networks,” _IEEE Trans. Aerosp. Electron. Syst._ , vol. 58, no. 5, pp. 3881–3894, Aug. 2022.
  * [12] M. Haenggi, _Stochastic Geometry for Wireless Networks_. Cambridge Univ. Press, 2012.
  * [13] M. Alzenad and H. Yanikomeroglu, “Coverage and rate analysis for vertical heterogeneous networks (VHetNets),” _IEEE Trans. on Wireless Commun._ , vol. 18, no. 12, pp. 5643–5657, Dec. 2019.
  * [14] I. Gordon _et al._ , “The HITRAN2016 molecular spectroscopic database,” _J. Quant. Spectrosc. Radiat. Transf._ , vol. 203, pp. 3–69, Dec. 2017.
  * [15] K. Belbase, Z. Zhang, H. Jiang, and C. Tellambura, “Coverage analysis of millimeter wave decode-and-forward networks with best relay selection,” _IEEE Access_ , vol. 6, pp. 22 670–22 683, 2018.

  *[SNR]: signal-to-noise ratio
  *[CCDF]: complementary cumulative distribution function
  *[PDF]: probability density functions


# Proprioceptive State Estimation of Legged Robots with Kinematic Chain
Modeling

Varun Agrawal1, Sylvain Bertrand2, Robert Griffin2,3, Frank Dellaert1
1Institute for Robotics and Intelligent Machines and School of Interactive
Computing, Georgia Institute of Technology, Atlanta, GA 30332 USA.
<EMAIL_ADDRESS>Institute for Human & Machine
Cognition, Pensacola, FL 32502 USA<EMAIL_ADDRESS>Systems and Robots, University of West Florida, Pensacola, FL 32503 USA.

###### Abstract

Legged robot locomotion is a challenging task due to a myriad of sub-problems,
such as the hybrid dynamics of foot contact and the effects of the desired
gait on the terrain. Accurate and efficient state estimation of the floating
base and the feet joints can help alleviate much of these issues by providing
feedback information to robot controllers. Current state estimation methods
are highly reliant on a conjunction of visual and inertial measurements to
provide real-time estimates, thus being handicapped in perceptually poor
environments. In this work, we show that by leveraging the kinematic chain
model of the robot via a factor graph formulation, we can perform state
estimation of the base and the leg joints using primarily proprioceptive
inertial data. We perform state estimation using a combination of
preintegrated IMU measurements, forward kinematic computations, and contact
detections in a factor-graph based framework, allowing our state estimate to
be constrained by the robot model. Experimental results in simulation and on
hardware show that our approach out-performs current proprioceptive state
estimation methods by 27% on average, while being generalizable to a variety
of legged robot platforms. We demonstrate our results both quantitatively and
qualitatively on a wide variety of trajectories.

##  I INTRODUCTION

Legged robots are considered advantageous in comparison to traditional wheeled
robots due to their ability to traverse uneven [1, 2] and unstructured
terrains [3]. Deployment of robots in the real-world is difficult since the
world is designed to be human-centric. Environments such as multi-level
buildings, and unstable or discontinuous ground [4, 5] are better tackled by
legged robots. Humanoids, in particular, are well suited for tasks which
require human levels of dexterity, manipulation, and characterization [6, 7,
8].

Many challenges with legged robot locomotion and stability can be addressed
with accurate and efficient state estimation. The robot’s high degree of
freedom from its articulated legs, the hybrid dynamics due to discrete foot-
contacts, the desired gait specifications, and managing foot slip [9] are some
of the key challenges which can be tackled by accurately estimating the state
of the robot’s floating base and the trajectory of foot poses. In particular,
high frequency state feedback controllers require even higher rate state
estimates, which is necessary to deal with uncertain events such as foot slip
and balance.

((a)) The DRC-era Atlas humanoid by Boston Dynamics in action.

((b)) The estimated trajectory of the Atlas walking through obstacles.

((c)) The estimated translation values versus the ground truth.

((d)) The estimated rotation values versus the ground truth.

Figure 1: A representative example of our state estimator working for the
Atlas robot. Experiments are conducted on the A1 robot in simulation and on
real hardware for the Anymal C robot by ANYbotics AG and the Atlas,
demonstrating the generality of our approach.

Current legged robot state estimators leverage exteroceptive sensors which
fail in poor environmental conditions and limit the frequency of controllers.
Exteroceptive sensors, such as cameras and LIDARs, are used in conjunction
with proprioceptive sensors such as Inertial Measurement Units (IMUs) and
joint encoders [10] [11] [12], and these methods are reliant on being able to
detect and process features in the surroundings. This requires well-lit,
opaque, and structured environments, which is not always the case, while also
being limited by the frame rate of visual sensors (e.g. cameras), which tend
to be an order of magnitude slower than that of IMUs, limiting the publishing
rate of state estimates, and subsequently the controller.

Proprioceptive state estimation has shown to be a viable solution [13, 14,
15], however it suffers from a fair number of problems. While having a
multitude of advantages, such as IMU-rate (e.g. 200 Hz) estimates for dynamic
control, operating in perception-denied environments [16], and freeing up
sensors (e.g. cameras) for other tasks such as mapping [11], they suffer from
long-term drift (particularly for filtering based approaches [13]), vary for
different legged robot platforms, and assume one noise model from the base to
the foot without taking into account complex joint interactions and link
deflections introduced over time.

We propose a framework to perform kinematic state estimation for legged robots
using factor graphs which tackles both issues of drift and generalizability.
Using proprioceptive measurements with the kinematic chain of the robot’s legs
in contact, available from the Unified Robot Description Format (URDF) file,
we can directly model the full kinematics of the robot at a time instance with
a factor graph, and perform probabilistic inference for the full state of the
robot. Moreover, the factor graphs allow for a generic sensor-fusion
mechanism, allowing the incorporation of additional factors and sensors such
as leg IMUs, contact sensors, and slip-detectors, depending on the situation,
as well as an extensible framework for asynchronous data processing.

Our method consistently produces better state estimates than the baseline, in
terms of Absolute Pose Error (APE) and Relative Pose Error (RPE), on a variety
of trajectories for each platform. We demonstrate our results in simulation
using the Unitree A1 robot [17, 2], and on physical hardware using the Anymal
C quadruped [18] and the Boston Dynamics’ DRC-era Atlas humanoid [9]. As a
comparative baseline, we use a factor graph based implementation of [19]. We
choose this baseline since it illustrates the benefit of our approach directly
in a fair comparison, while being orthogonal to more recent advances.

##  II RELATED WORK

Use of inertial measurement units (IMUs) for proprioceptive state estimation
has a rich history [20]. The initial application was for the control of
intercontinental ballistic missiles [21]. Being mission-critical, they are
highly sophisticated (tactical or navigation grade) but also very expensive
[21, Table 1], making them typically unsuitable for general robotics
applications.

More recently, IMUs have seen widespread adoption in the areas of visual-
inertial odometry and robot localization [22, 23, 24]. Initial work focused on
sensor fusion based on Kalman filtering [25, 26, 27, 28]. Recent smoothing
methods, such as [29], leverage preintegrated IMU measurements with camera-
based visual measurements to perform robot localization [30] and state
estimation [31]. While these are synergistic with goals such as mapping, they
suffer from issues such as lighting dependency, potential lack of visual
features, and increased computational cost, causing up to two orders of
magnitude slower processing rates.

For legged robots, initial work with proprioceptive state estimation was
performed using Kalman filtering and the contact points of the robot with its
environment [32, 19, 13]. These foundational works focus on using the forward
kinematics of the legs to define constraints on the position and velocity of
the contact points and thus constrain the biases in the IMU measurements. The
major drawback is that the yaw angle (around gravity) and the absolute
position of the robot become unobservable, leading to drift errors over long
durations. Related to our work, [33] perform filtering based state estimation
for humanoids specifically, using ideas from [19] while leveraging the flat-
foot constraint of these robots to tackle yaw angle unobservability.

Recently, there has been an increased interest in smoothing based methods for
legged robot state estimation. [34] demonstrate excellent results for state
estimation by fusing proprioceptive and exteroceptive sensor measurements in a
factor-graph based framework for a quadruped. However, for leg odometry, they
rely on interpolation over the states from a Two-State Implicit Filter [35]
rather than jointly estimate it. Closely related to our work, [36] leverage a
factor-graph smoothing framework to jointly optimize over IMU measurements,
leg kinematics and contact measurements preintegrated over each stance phase
of biped locomotion. A follow up work [12] shows the promising potential for
hybrid system modeling of contacts inherent to bipedal walking. By adding
these hybrid preintegrated contact factors with visual features from a stereo
camera setup, they are able to demonstrate state estimation over a variety of
conditions and terrains.

##  III PRELIMINARIES

In this section, we present some preliminaries on factor graphs, the robot
kinematic model, and IMU preintegration in the context of smoothing, all of
which serve as the basis for the discussion of our work in subsequent
sections.

###  III-A Factor Graph Based Smoothing

A factor graph is a probabilistic graph model which can be used to express the
estimation problem in a simple and elegant manner. It is a bipartite graph
represented as $\textit{G}=\\{\phi,\Theta,\varepsilon\\}$, where $\phi$ are
the factor nodes, $\Theta$ are the variable nodes to be estimated, and
$\varepsilon$ are the edges between them.

The factor graph G represents the probability density over the variables to be
estimated $\Theta$ given the measurement values $Z$, and provides the
factorization as

$P(\Theta|Z)=\prod_{i=1}^{|\textit{G}|}\phi_{i}(\Theta_{i}),$ (1)

with each factor $\phi_{i}(\Theta_{i})=P(\Theta_{i}|Z_{i})$ specifying the
relationship between its connected subset of variables
$\Theta_{i}\subseteq\Theta$.

Given the factor graph G, the objective is to compute the maximum a posteriori
estimate of variables $\Theta$, given as

$\Theta^{*}=\arg\max_{\Theta}\prod_{i=1}^{|\textit{G}|}\phi_{i}(\Theta_{i})=\arg\min_{\Theta}-\sum_{i=1}^{|\textit{G}|}\text{log}\phi_{i}(\Theta_{i}).$
(2)

###  III-B IMU Preintegration

In this section, we provide a brief overview of IMU preintegration and how it
is used in the factor graph framework. For more details, the reader is
referred to the exposition in [29].

It is important to handle an IMU’s high measurement frequency so as to not
overwhelm the state estimator when performing smoothing. This is done by
forward integrating the IMU’s measurements (angular velocity $\omega$ and
linear acceleration $a$) across successive $\Delta t$ between two specific
timestamps to obtain a summarized compound measurement for the body $B$ in an
arbitrary world frame $W$. To prevent the forward integration being performed
repeatedly at each new linearization point, a relative motion model is
computed between the states at time $t_{i}$ and time $t_{j}$

$\displaystyle\begin{split}&\Delta\mathrm{R}_{ij}\triangleq\prod_{k=i}^{j-1}\text{Exp}\left(\left(\omega_{k}-\mathrm{b}^{g}_{k}-\eta^{gd}_{k}\right)\Delta
t\right),\\\
&\Delta\mathrm{v}_{ij}\triangleq\sum_{k=i}^{j-1}\Delta\mathrm{R}_{ik}\left(a_{k}-\mathrm{b}^{a}_{k}-\eta^{ad}_{k}\right)\Delta
t,\\\
&\Delta\mathrm{p}_{ij}\triangleq\sum_{k=i}^{j-1}\left[\Delta\mathrm{v}_{ik}\Delta
t+\frac{1}{2}\Delta\mathrm{R}_{ik}(a_{k}-\mathrm{b}^{a}_{k}-\eta_{k}^{ad}\Delta
t^{2})\right],\end{split}$ (3)

where $\Delta\mathrm{R}_{ij}$, $\Delta\mathrm{v}_{ij}$ and
$\Delta\mathrm{p}_{ij}$ are the rotation, velocity and translation differences
between times $t_{i}$ and $t_{j}$ denoted by $\Delta t$, $ad$ is the
accelerometer standard deviation and $gd$ is the gyroscope standard deviation.

This summarized relative state is used as a motion model constraint between
the previous and current state of the body. By isolating the noise terms, we
get

$\displaystyle\begin{split}&\Delta\mathrm{R}_{ij}=\mathrm{R}_{i}^{T}\mathrm{R}_{j}\text{Exp}(\delta\phi_{ij}),\\\
&\Delta\mathrm{v}_{ij}=\mathrm{R}_{i}^{T}(\mathrm{v}_{j}-\mathrm{v}_{i}-\text{g}\Delta
t_{ij})+\delta\mathrm{v}_{ij},\\\
&\Delta\mathrm{p}_{ij}=\mathrm{R}_{i}^{T}(\mathrm{p}_{j}-\mathrm{p}_{i}-\mathrm{v}_{i}\Delta
t_{ij}-\frac{1}{2}\text{g}\Delta
t_{ij}^{2})+\delta\mathrm{p}_{ij}.\end{split}$ (4)

with $g$ being the value of Earth’s gravity, and $\delta\phi_{ij}$,
$\delta\mathrm{v}_{ij}$ and $\delta\mathrm{p}_{ij}$ being the rotation,
velocity and translation increments between times $t_{i}$ and $t_{j}$ due to
the noise. For details on noise and covariance propagation, please refer to
[29].

##  IV APPROACH

In this section, the state estimation problem is formulated and we describe
the overall approach modeled by the factor graph. We assume the IMU is rigidly
attached to the base link, and that the robot either has force/tactile sensors
on its feet or has some contact detection mechanism [37] to give us foot
contact information.

###  IV-A Problem Formulation

The state variables are denoted by the pose of the base link $P^{W}_{B}$ in
some arbitrary world coordinate frame $W$, and its linear velocity $V_{B}$.
Each foot link has its pose at contact in the world frame denoted by
$P^{W}_{L_{i}}$, where $i\in\\{N\\}$ is the index of the corresponding leg for
$N$ legs. We denote the bias of the accelerometer and the gyroscope as
$\mathrm{b}^{a}$ and $\mathrm{b}^{g}$ respectively. This gives the robot state
at time $t_{i}$ as

$x_{i}\dot{=}[P^{W}_{Bi},V_{B},\mathrm{b}^{g},\mathrm{b}^{a}]\in\textbf{SE}(3)\times\mathbb{R}^{9}.$
(5)

Between two time steps $i$ and $j$, we denote the IMU measurements as
$(\omega_{i,j},a_{i,j})$, the joint measurements as $J_{i,j}$ and the contact
states for all legs as $C_{i,j}$. The set of all measurements up to time-step
$k$ is denoted by

$Z^{k}\dot{=}\\{\omega^{k},a^{k},J^{k},C^{k}\\}.$ (6)

Given the measurements from the IMU, joint encoders, and contact sensors, we
frame our problem as probabilistic inference. By assuming the measurements are
conditionally independent and the noise on the sensors are zero mean additive
Gaussians, the posterior probability of the state estimation problem can be
stated as

$\displaystyle\phi(X^{K})\propto p(X^{K}|Z^{k})$
$\displaystyle\dot{=}\phi(x_{0})\prod_{k=1}^{K-1}\phi^{\text{IMU}}(x_{k},x_{k-1})\prod_{j=1}^{L}\phi^{FK}(x_{k},l_{k,j})\phi^{CP}(l_{k,j},c_{i})$

where $\phi^{IMU}$ are the preintegrated IMU factors, $\phi^{FK}$ are the
Forward Kinematics factors, and $\phi^{CP}$ are the factors equating the leg
end-effector to the same contact point through a stance phase.

Figure 2: The factor graphs for the baseline estimator and our proposed state
estimator. We propose modeling each joint in the leg’s kinematic chain to
apply a noise model, removing the need for noise assumptions on the forward
kinematics between the base and the end link. Here, in the case of the A1, B
represents the base link, T the thigh link, U the upper leg link, and L the
lower leg link. The variable for the contact point is denoted by CP.

###  IV-B Preintegrated IMU Measurements

In the same vein as [19, 13], we use the measurement data from the IMU located
at the base of the robot. We perform preintegration of the IMU measurements
[29] and add a preintegrated IMU factor at the specified estimation rate (50
Hz for quadrupeds, 100 Hz for the humanoid). The estimation rate for factor
addition can be updated depending on the use case, allowing us to recover
estimates at a rate suitable for use in trajectory controllers. A smoothing
based approach also helps avoid drift due to accumulated linearization errors,
commonly seen in filtering based approaches.

###  IV-C Kinodynamic Characterization with Factor Graphs

The robot kinematic state can be modeled as a factor graph as the links and
joints follow a similar bipartite structure [38]. A robot’s kinematic model
can be specified as $\textit{R}=\\{J,L,\mathfrak{E}\\}$ which is also a
bipartite graph between joints $J$ and links $L$, with the edges
$\mathfrak{E}$ specifying the connections. Thus, for each joint we can specify
a joint measurement factor $\phi^{FK}_{j}$ which takes the joint angle as a
measurement and provides the relationship between the connecting links
$\\{L_{j},L_{j+1}i\\}$.

An additional, non-obvious, benefit of the factor graph framework is the
potential for extensibility. One can add additional factors to incorporate
prior knowledge as well as encode additional properties of the robot such as
its dynamics (e.g. inertial matrix, friction cone). This allows for
potentially also estimating the centroidal dynamics of the robot [39] which we
leave for future work.

###  IV-D Forward Kinematics at Contact

At each point a leg is in contact, we add the sub-graph to model the kinematic
chain from the base of the robot to the link in contact. This allows us to
probabilistically model the interactions between the links’ noise models,
rather than assume gaussian noise over the forward kinematics. In contrast to
[13, 36] which uses a single forward kinematics constraint between the base
link and the foot link, we add multiple pose factors which enforce a
transformation between each joint’s parent and child links. Each factor takes
as measurement the angle of the joint as recorded by a joint encoder or joint-
level IMU. By using a minimal representation of the joint’s screw axis as
described in [40], we are able to ensure an efficient representation despite
the addition of more factors. Fig 2 provides a visualization.

For $J^{f}$ joints in the kinematic chain for each leg $f$ in the set of legs
$F$, each joint encoder measurement $z_{k,j}$, and contact $c_{k,f}$ at
timestep $k$, we get the following factorization of the posterior
$p(X_{k},L_{k,1:F}|z_{k,1:J},c_{k,1:F})$ on the base state $X_{k}$ (equal to
the link pose $L^{1}_{k}$) and all leg links $L_{k,1:F}$:

$\displaystyle p(X_{k},L_{k,1:F}|z_{k,1:J},c_{k,1:F})\propto$ (7)
$\displaystyle\prod_{f\in
F}\phi^{FK}(X_{k},L^{2}_{k,f}|z_{k,1},c_{k,f})\prod^{J^{f}}_{j=2}\phi^{FK}(L^{j}_{k,f},L^{j+1}_{k,f}|z_{k,j},c_{k,f})$

###  IV-E Contact Points as Landmarks

In order to constrain the foot positions and estimate the IMU biases, we also
estimate the contact points as being consistent across the leg’s stance phase.
By viewing the state estimation problem as a SLAM problem, we consider the
contact points as landmarks to be estimated and use techniques from SLAM (e.g.
data association) directly.

The Forward Kinematics ($\mathrm{FK}$) provides a mapping between the base
link of the robot and the foot contact point, given that the foot is in
contact. We specify a foot is in contact with point $m\in\\{M\\}$ at time-step
$k$ using the factor $\phi^{CP}(X_{k},C_{im})$ which minimizes the distance
between the contact point $m$ in the environment and the foot link $i\in F$

$\displaystyle\phi^{CP}(X_{k},C_{im})=\mathrm{g}^{CP}(X_{k},C_{im})=\|\mathrm{FK}(X_{k},L_{i})-C_{im}\|$
(8)

Depending on the type of foot, either Point or Flat, we can constrain just the
foot translation or the foot pose [33].

###  IV-F Maximum A Posteriori Estimation

Given a linearization point $X$, we linearize our factor graph and optimize
for the negative log-likelihood to give us the best estimates of the state
trajectory using the Levenberg-Marquardt optimization method

$\displaystyle
X^{*}=\arg\min_{X}\bigg{(}\|o(X_{0})\|^{2}_{\Sigma_{0}}+\sum_{k=1}^{K-1}\|\mathrm{f}^{\mathrm{IMU}}(X_{k},X_{k-1})\|^{2}_{\Sigma_{\mathrm{IMU}}}$
(9)
$\displaystyle+\sum_{k=0}^{K-1}\sum_{f}^{F}\|\mathrm{h}^{\mathrm{FK}}(X_{k},L_{k,f}^{2})\|^{2}_{\Sigma_{\mathrm{FK}}}\sum_{i=2}^{N}\|\mathrm{h}^{\mathrm{FK}}(L^{i}_{k,f},L^{i+1}_{k,f})\|^{2}_{\Sigma_{\mathrm{FK}}}$
$\displaystyle+\sum_{k=0}^{K-1}\sum_{i}^{F}\sum_{m}^{M}\|\mathrm{g}^{\mathrm{CP}}(X_{k},C_{im})\|^{2}_{\Sigma_{\mathrm{CP}}}\bigg{)}$

where $o(X_{0})$ is a prior on the first state, $\mathrm{f}^{\mathrm{IMU}}$,
$\mathrm{h}^{\mathrm{FK}}$ and $\mathrm{g}^{\mathrm{CP}}$ are the linearized
binary preintegrated IMU error, forward kinematics error for a parent-child
link pair, and the contact error for foot $i$, respectively.
$\Sigma_{0},\Sigma_{\mathrm{IMU}},\Sigma_{\mathrm{FK}},\Sigma_{\mathrm{CP}}$
are the respective noise model covariances.

We compute a good linearization point using the preintegrated IMU measurements
and the joint angle measurements, giving us initial values for the state and
the contact points.

##  V Experimental Results

To evaluate our approach, we use common metrics used for SLAM[41, 42] such as
Absolute Pose Error and Relative Pose Error for investigating the global and
local consistency of the trajectories respectively. We make use of the evo
software package [43] to provide us with convenient functions for the metrics.
The Absolute Pose Error (APE) is defined as the absolute relative pose between
two poses, the reference pose $P_{i}$ and the estimated pose $\hat{P}_{i}$ at
timestamp $i$

$APE_{i}=\|(P_{i}^{-1}\hat{P}_{i})-I_{4\times 4}\|_{F}$ (10)

The Relative Pose Error (RPE) compares the relative poses along the estimated
and reference trajectories and is computed over $1$ second

$RPE_{i,j}=\|[(P_{i}^{-1}P_{j})^{-1}(\hat{P}_{i}^{-1}\hat{P}_{j})]-I_{4\times
4}\|_{F}$ (11)

To compute the metrics, we first have to align the trajectories, which we do
with a novel scheme for legged robots. Trajectory alignment ensures that the
error being computed is not dominated by any particular unobservable gauge
freedom, e.g. scale in traditional SLAM systems. In our case, we have 4
degrees of gauge freedom, namely the translation (contributing 3 degrees of
freedom) and the yaw $\gamma$, all of which are directly unobservable for
legged robots.

The trajectory alignment is done in two steps, the $(x,y,\gamma)$ and the $z$
axis. Aligning the $(x,y,\gamma)$ is a straightforward $SE(2)$ transformation
calculated by computing the best alignment between the $(x,y)$ points of the
trajectories. For aligning the $z$ axis, we subtract the mean of the $z$
values of the trajectory translations from each trajectory independently.

TABLE I: Experimental Results on Quadrupeds in both simulation (A1) and hardware (Anymal C). We report the Root Mean Square Error (RMSE) for Absolute Pose Error (APE) and Relative Pose Error (RPE) over 1 second. | Baseline Estimator | Our Approach
---|---|---
Trajectory | APE | RPE | APE | RPE
A1 Straight | 0.207 | 0.0063 | 0.205 | 0.0063
A1 Diagonal | 0.136 | 0.0064 | 0.135 | 0.0064
A1 Turn | 0.190 | 0.0053 | 0.181 | 0.0049
A1 Zig-Zag | 0.109 | 0.0062 | 0.108 | 0.0063
Anymal C | 0.307 | 0.0115 | 0.173 | 0.0046
TABLE II: Experimental Results on a DRC-era Atlas humanoid over 5 different test runs. Similar to Table I, we report the Root Mean Square Error (RMSE) for Absolute Pose Error (APE) and Relative Pose Error (RPE) over 1 second for each trajectory. | Baseline Estimator | Our Approach
---|---|---
Trajectory | APE | RPE | APE | RPE
Atlas Run 1 | 0.076 | 0.0201 | 0.067 | 0.0182
Atlas Run 2 | 0.192 | 0.0339 | 0.141 | 0.0277
Atlas Run 3 | 0.198 | 0.0313 | 0.098 | 0.0201
Atlas Run 4 | 0.154 | 0.0267 | 0.114 | 0.0189
Atlas Run 5 | 0.292 | 0.0324 | 0.271 | 0.0260

We evaluate our estimator in simulation to show correctness and establish the
strong baseline. We use the pybullet based simulation environment provided by
[17]. We generate four different trajectories in simulation which we refer to
as straight, diagonal, turn, and zig-zag. In the simulation, along with the
ground truth state values, we collect the joint angles for all the joints in
the robot and the contact state of each leg. We add realistic noise to the
simulated measurements, i.e., from the IMU and the joint angles. We execute
both the baseline estimator and our estimator on the different trajectories.

The baseline estimator used is a factor-graph based implementation of the
proprioceptive state estimator described by [13]. We choose this baseline
since it is an estimator that is widely deployed in many real-world robots
[44, 14, 15]. By implementing a factor graph based version of the original
filter-based estimator, we bring forward the benefits of the original filter
while also allowing for a fair comparison between the baseline and our
proposed framework.

Figure 3: The A1 walking a zig-zag trajectory in simulation. Qualitative
comparisons of (a) the trajectories, (b) the translations, and (c) the
rotations, verifying the performance of our state estimator in ideal
conditions.

Finally, we demonstrate the efficacy of our estimator on robot hardware. We
run our estimator on a sample trajectory of the ANYbotics Anymal C robot and a
Boston Dynamics Atlas robot. We track the trajectory of the Anymal C base
using a Vicon optical tracking system, akin to [10]. For the Atlas, we use the
estimates provided by the well-tuned built-in state estimator as a pseudo-
ground-truth trajectory.

Quantitative results for the A1 simulation and the Anymal C in Table I and for
the Atlas in Table II, show our state estimator’s superior performance over a
variety of trajectories for both platforms. As seen from the reported results,
our state estimator outperforms the baseline estimator in APE and RPE metrics
by 27.62% and 28.75% on average respectively. The RPE improvement illustrates
improved drift handling over 1 second intervals, and in particular, for 5
different trajectories of the Atlas humanoid, we show improved state estimates
consistently, demonstrating both the accuracy and the capability of our
approach. Sample trajectories are displayed in Figures 3, 4, 5.

Figure 4: State estimation for the Anymal C robot walking a non-trivial
trajectory. The ground truth trajectory was gathered using motion capture
technology in a lab environment.

Figure 5: State estimation results for the Atlas robot. Both rotation and
translation estimates are equivalent or outperform the baseline estimator.

##  VI CONCLUSION

In this work, we presented a factor-graph based framework for proprioceptive
legged robot state estimation. By viewing the robot’s kinematic structure as a
bipartite graph, we are able to model it using a factor graph together with
the IMU preintegration factors. This in-turn allows for constraining the IMU
biases for accurate and efficient proprioceptive state estimation. We
demonstrate the efficacy and generalizability of our approach on both
humanoids and quadrupeds, showing improved performance when compared to a
strong baseline in both simulation and robot hardware.

Our framework is extensible and allows for asynchronous data handling. In
future work, we hope to examine the effects of new factors based on force-
torque measurements and constraints based on slip-estimation. Our hope is that
this work serves as a stepping-stone towards more widespread deployment of
legged robots for various applications.

## ACKNOWLEDGMENTS

We thank Russell Buchanan and Bhavyansh Mishra for help with the Anymal C
robot data and the Atlas robot data, respectively. We would also like to thank
Jerry Pratt, Maurice Fallon and Fan Jiang for helpful discussions and advice.

## References

  * [1] Ruben Grandia, Farbod Farshidian, René Ranftl and Marco Hutter  “Feedback MPC for Torque-Controlled Legged Robots”  In _2019 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2019, pp. 4730–4737  DOI: 10.1109/IROS40897.2019.8968251
  * [2] Ashish Kumar, Zipeng Fu, Deepak Pathak and Jitendra Malik  “RMA: Rapid Motor Adaptation for Legged Robots”, 2021  arXiv:2107.04034 [cs.LG]
  * [3] Hendrik Kolvenbach et al.  “Towards autonomous inspection of concrete deterioration in sewers with legged robots”  In _Journal of Field Robotics_ 37.8, 2020, pp. 1314–1327  DOI: https://doi.org/10.1002/rob.21964
  * [4] V. Lebastard, Y. Aoustin and F. Plestan  “Estimation of Absolute Orientation for a Bipedal Robot: Experimental Results”  In _IEEE Transactions on Robotics_ 27.1, 2011, pp. 170–174  DOI: 10.1109/TRO.2010.2094410
  * [5] Claudio Semini et al.  “HyQ - a dynamic locomotion research platform”, 2011
  * [6] C. Breazeal  “Emotion and sociable humanoid robots”  In _Int. J. Hum. Comput. Stud._ 59, 2003, pp. 119–155
  * [7] K. Hirai, M. Hirose, Y. Haikawa and T. Takenaka  “The development of Honda humanoid robot”  In _Proceedings. 1998 IEEE International Conference on Robotics and Automation (Cat. No.98CH36146)_ 2, 1998, pp. 1321–1326 vol.2
  * [8] Gabe Nelson et al.  “PETMAN: A Humanoid Robot for Testing Chemical Protective Clothing”  In _Journal of the Robotics Society of Japan_ 30.4, 2012, pp. 372–377  DOI: 10.7210/jrsj.30.372
  * [9] Matthew Johnson et al.  “Team IHMC’s Lessons Learned from the DARPA Robotics Challenge Trials”  In _J. Field Robot._ 32.2 GBR: John WileySons Ltd., 2015, pp. 192–208  DOI: 10.1002/rob.21571
  * [10] David Wisth, Marco Camurri and Maurice Fallon  “VILENS: Visual, Inertial, Lidar, and Leg Odometry for All-Terrain Legged Robots”, 2021  arXiv:2107.07243 [cs.RO]
  * [11] M Camurri, M Ramezani, S Nobili and M Fallon  “Pronto: A multi-sensor state estimator for legged robots in real-world scenarios”  In _Frontiers in Robotics and AI_ 7 Frontiers in RoboticsAI, 2020
  * [12] Ross Hartley et al.  “Hybrid Contact Preintegration for Visual-Inertial-Contact State Estimation Using Factor Graphs”  In _Proceedings of the IEEE/RSJ International Conference on Intelligent Robots and Systems_ , 2018, pp. 3783–3790  URL: https://arxiv.org/pdf/1803.07531.pdf
  * [13] M. Bloesch et al.  “State Estimation for Legged Robots: Consistent Fusion of Leg Kinematics and IMU”  In _Robotics: Science and Systems VIII_ IEEE, 2013, pp. 17–24
  * [14] Gerardo Bledt et al.  “MIT Cheetah 3: Design and Control of a Robust, Dynamic Quadruped Robot”  In _2018 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2018, pp. 2245–2252  DOI: 10.1109/IROS.2018.8593885
  * [15] M. Fallon A. Valenzuela H. Dai F. Permenter T.ΩKoolen P. Marion S. Kuindersma and R. Tedrake  “Optimization-based locomotion planning, estimation, and control design for the atlas humanoid robot”  In _Autonomous robots_ 40.3 Springer, 2016, pp. 429–455
  * [16] Marco Tranzatto et al.  “CERBERUS: Autonomous Legged and Aerial Robotic Exploration in the Tunnel and Urban Circuits of the DARPA Subterranean Challenge”  arXiv, 2022  DOI: 10.48550/ARXIV.2201.07067
  * [17] Xue Bin Peng et al.  “Learning Agile Robotic Locomotion Skills by Imitating Animals”  In _ArXiv_ abs/2004.00784, 2020
  * [18] David Wisth, Marco Camurri and Maurice Fallon  “Preintegrated Velocity Bias Estimation to Overcome Contact Nonlinearities in Legged Robot Odometry”  In _IEEE Intl. Conf. on Robotics and Automation (ICRA)_ , 2020
  * [19] M. Bloesch et al.  “State estimation for legged robots on unstable and slippery terrain”  In _2013 IEEE/RSJ International Conference on Intelligent Robots and Systems_ , 2013, pp. 6058–6064  DOI: 10.1109/IROS.2013.6697236
  * [20] Robert D. Braun et al.  “Advances in Inertial Guidance Technology for Aerospace Systems”  In _AIAA Guidance, Navigation, and Control (GNC) Conference_ , 2013  DOI: 10.2514/6.2013-5123
  * [21] Scott M. Bezick, Alan J. Pue and Charles M. Patzelt  “Inertial Navigation for Guided Missile Systems”  In _Johns Hopkins Apl Technical Digest_ 28, 2010, pp. 331–342
  * [22] S. Leutenegger et al.  “Keyframe-based visual-inertial odometry using nonlinear optimization”  In _Intl. J. of Robotics Research_ 34.3, 2015, pp. 314–334
  * [23] Raul Mur-Artal and Juan D. Tardos  “Visual-Inertial Monocular SLAM with Map Reuse” arXiv:1610.05949, 2016
  * [24] T. Qin, P. Li and S. Shen  “VINS-Mono: A Robust and Versatile Monocular Visual-Inertial State Estimator”  In _IEEE Transactions on Robotics_ 34.4, 2018, pp. 1004–1020  DOI: 10.1109/TRO.2018.2853729
  * [25] M. Bloesch, S. Omari, M. Hutter and R. Siegwart  “Robust visual inertial odometry using a direct EKF-based approach”  In _2015 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2015, pp. 298–304  DOI: 10.1109/IROS.2015.7353389
  * [26] A.I. Mourikis and S.I. Roumeliotis  “A Multi-State Constraint Kalman Filter for Vision-aided Inertial Navigation”  In _IEEE Intl. Conf. on Robotics and Automation (ICRA)_ , 2007, pp. 3565–3572
  * [27] V. Usenko, J. Engel, J. Stückler and D. Cremers  “Direct Visual-Inertial Odometry with Stereo Cameras”  In _IEEE Intl. Conf. on Robotics and Automation (ICRA)_ , 2016, pp. 1885–1892
  * [28] Stephan Weiss et al.  “Real-time onboard visual-inertial state estimation and self-calibration of mavs in unknown environments”  In _IEEE Intl. Conf. on Robotics and Automation (ICRA)_ , 2012, pp. 957–964  IEEE
  * [29] C. Forster, L. Carlone, F. Dellaert and D. Scaramuzza  “On-Manifold Preintegration Theory for Fast and Accurate Visual-Inertial Navigation”  In _IEEE Trans. Robotics_ 33.1, 2017, pp. 1–21
  * [30] Kevin Eckenhoff, Patrick Geneva and Guoquan Huang  “Closed-form preintegration methods for graph-based visual–inertial navigation”  In _The International Journal of Robotics Research_ 38.5, 2019, pp. 563–586  DOI: 10.1177/0278364919835021
  * [31] Barza Nisar, Philipp Foehn, Davide Falanga and D. Scaramuzza  “VIMO: Simultaneous Visual Inertial Model-Based Odometry and Force Estimation”  In _IEEE Robotics and Automation Letters_ 4, 2019, pp. 2785–2792
  * [32] M. Bloesch and M. Hutter  “Technical Implementations of the Sense of Balance”  In _Humanoid Robotics: A Reference_ Dordrecht: Springer Netherlands, 2016, pp. 1–29  DOI: 10.1007/978-94-007-7194-9˙69-2
  * [33] N. Rotella, M. Bloesch, L. Righetti and S. Schaal  “State estimation for a humanoid robot”  In _2014 IEEE/RSJ International Conference on Intelligent Robots and Systems_ , 2014, pp. 952–958  DOI: 10.1109/IROS.2014.6942674
  * [34] David Wisth, Marco Camurri and Maurice Fallon  “Robust Legged Robot State Estimation Using Factor Graph Optimization”  In _IEEE Robotics and Automation Letters_ , 2019
  * [35] Michael Bloesch et al.  “The Two-State Implicit Filter - Recursive Estimation for Mobile Robots”  In _IEEE Robotics and Automation Letters_ 3.1, 2017, pp. 573–580  DOI: 10.1109/LRA.2017.2776340
  * [36] Ross Hartley et al.  “Legged Robot State-Estimation Through Combined Forward Kinematic and Preintegrated Contact Factors”  In _Proceedings of the IEEE International Conference on Robotics and Automation_ , 2018, pp. 4422–4429  URL: https://arxiv.org/pdf/1712.05873.pdf
  * [37] Fabian Jenelten et al.  “Dynamic Locomotion on Slippery Ground”  In _IEEE Robotics and Automation Letters_ 4.4, 2019, pp. 4170–4176  DOI: 10.1109/LRA.2019.2931284
  * [38] Mandy Xie and Frank Dellaert  “Batch and Incremental Kinodynamic Motion Planning using Dynamic Factor Graphs”  In _arXiv preprint arXiv:2005.12514_ , 2020
  * [39] Hongkai Dai, Andrés Valenzuela and Russ Tedrake  “Whole-body motion planning with centroidal dynamics and full kinematics”  In _2014 IEEE-RAS International Conference on Humanoid Robots_ , 2014, pp. 295–302  DOI: 10.1109/HUMANOIDS.2014.7041375
  * [40] K. Lynch and F. Park  “Modern Robotics”  Cambridge University Press, 2017  URL: https://dl.acm.org/doi/book/10.5555/3165183
  * [41] Jürgen Sturm et al.  “A benchmark for the evaluation of RGB-D SLAM systems”  In _2012 IEEE/RSJ International Conference on Intelligent Robots and Systems_ , 2012, pp. 573–580  DOI: 10.1109/IROS.2012.6385773
  * [42] R. Kümmerle et al.  “On Measuring the Accuracy of SLAM Algorithms”  In _Autonomous Robots_ 27.4, 2009, pp. 387–407  DOI: 10.1007/s10514-009-9155-6
  * [43] Michael Grupp  “evo: Python package for the evaluation of odometry and SLAM.”, https://github.com/MichaelGrupp/evo, 2017
  * [44] Marco Hutter et al.  “StarlETH: a Compliant Quadrupedal Robot for Fast, Efficient, and Versatile Locomotion”, 2012, pp. 483–490  DOI: 10.1142/9789814415958˙0062