# Accurate outdoor ground truth based on total stations

Maxime Vaidis, Philippe Giguère, François Pomerleau, Vladimír Kubelka
Northern Robotics Laboratory, Université Laval, Québec City, Québec, Canada
$\\{$maxime.vaidis, philippe.giguere, francois.pomerleau,
<EMAIL_ADDRESS>
###### Abstract

In robotics, accurate ground-truth position fostered the development of
mapping and localization algorithms through the creation of cornerstone
datasets. In outdoor environments and over long distances, total stations are
the most accurate and precise measurement instruments for this purpose. Most
total station-based systems in the literature are limited to three Degrees Of
Freedoms, due to the use of a single-prism tracking approach. In this paper,
we present preliminary work on measuring a full pose of a vehicle, bringing
the referencing system to six DOFs. Three total stations are used to track in
real time three prisms attached to a target platform. We describe the
structure of the referencing system and the protocol for acquiring the ground
truth with this system. We evaluated its precision in a variety of different
outdoor environments, ranging from open-sky to forest trails, and compare this
system with another popular source of reference position, the Real Time
Kinematics (RTK) positioning solution. Results show that our approach is the
most precise, reaching an average positional error of 10 mm and 0.6 deg. This
difference in performance was particularly stark in environments where Global
Navigation Satellite System (GNSS) signals can be weaker due to overreaching
vegetation.

###### Index Terms:

total station, motion tracking, ground truth, mobile robot, outdoor

##  I Introduction

Access to precise and accurate ground-truth pose reference during experiments
is vital in mobile robotics. For instance, knowing the pose of a robot, that
is its position and orientation, facilitates development and testing of
localization, mapping and control algorithms. For indoor experiments, motion
capture such as the Vicon [1, 2] or Optitrack [3] systems have become the _de
facto_ standard. They are based on tracking reflective or active markers by
infrared cameras and offer an accuracy of sub-millimetre precision and high
sampling rates. However, the area that they can cover is limited by the number
of cameras available. Moreover, they struggle in direct sunlight. Another
high-fidelity source of reference pose is laser trackers based on laser
interferometry, such as the one used by [4]. In exchange for the limitation of
tens of meters in range, they offer position measurements accurate down to
microns. The possibilities of pose reference are more limited outdoors. The
approaches most prominent in the literature use GNSS [5], total stations [6],
and qualitative comparison of trajectories to maps [7]. Besides, these outdoor
referencing systems do not provide direct orientation measurements.

Figure 1: Example of an application where RTK positioning faces limitations. A
robot navigating in a subarctic forest with three prismatic reflectors is
shown. It is being tracked by three total stations.

When it comes to gathering datasets in forests, such as the one depicted in
Figure 1, accuracy of RTK positioning suffers from signal attenuation by
vegetation [8]. Consequently, the need for accurate and precise pose reference
immune to foliage interference and signal reflection has motivated the
development of a system based on multiple total stations appropriate for such
task. Here, we make the distinction between a theodolite, which is a precision
optical measurement instrument used for measuring horizontal and vertical
angles, and a total station, which is a theodolite with distance-measuring
equipment.

In this article, we propose a protocol for a novel referencing system relying
on multiple total stations to allow the measurement of the six DOFs of a
vehicle. More precisely, our contributions are :

  1. 1.

the calibration procedure of multiple total stations for a vehicle continuous
tracking;

  2. 2.

the time synchronization and long-range data transmission protocol for online
pose referencing;

  3. 3.

a precision analysis of our tracking procedure; and

  4. 4.

a comparison against a RTK positioning system.

##  II Related work

In the robotics literature, we find a large variety of position-referencing
systems based on total stations. They generally differ in the types of
environment, in the number of referenced DOF, or in the target robotic systems
they are used for.

Regarding the robotic platforms, there are almost no limitations; the only one
being how many prisms the robot can carry. Researchers used total stations to
track skid steered robots [9, 10], a tethered wheeled robot [11], planetary
rover [12, 13], unmanned surface vessel [14], and walking robots [4, 15]. In
our application, we use a large skid-steered robot suitable for winter in
subarctic forest. However, the presented referencing system is generic and
applicable to any robot that can carry three prisms distanced by more than
$80\text{\,}\mathrm{cm}$.

Tracking a robot in a forest using a total station was previously done by [6].
They referenced the position of a wheeled tethered robot by a total station on
steep rocky slopes in a forested area to validate the accuracy of localization
and mapping algorithms. However, total stations are also resilient to harsh
environment, as they were used in deserts to track planetary rover [13].
Because of their long-range measurements, application on water bodies are
possible. For example, [14] used a total station installed on the shoreline to
track the position of a robotic catamaran. Total stations were also
demonstrated in urban areas [16] as well as indoors [17]. Given their
versatility in different environments, a referencing system relying on total
stations is thus relevant for multiple applications. The main limitation is
the requirement of the line-of-sight between the prisms and the total
stations.

An important design choice is the number of DOF the referencing system should
track. A total station naturally provides the position of a tracked target,
thus three DOFs. For example, in [18] they mounted a prismatic reflector, or
simply a prism, onto a race car to verify the output of its mapping and
localization system. Limited by their maximum payload capacity, system
tracking Unmanned Aerial Vehicles are limited to one prism [19, 20, 21]. An
interesting solution was presented by [22], which involves an actively-
controlled pan-tilt reflector on a climbing inspection robot. By always
pointing the reflector towards the total station, two DOFs were added to the
typical position measurement. Finally, with at least three prisms on a robot,
the full pose of a robot can be recovered using six DOFs. In the work of [23],
a lidar dataset was collected by a sensor rig with ground truth information
about the pose. A single total station was used to manually measure three
reflectors while the rig remained static. The same principle was used by [13]
for a planetary rover. We extend these works by proposing a system using
multiple total stations to continuously track a vehicle, while it is moving,
to recover its full pose.

##  III Referencing system design

###  III-A Referencing system overview

Total stations are better known for topographic surveys, but modern total
stations offer the ability to track a target over one kilometre with a range
precision within millimetres. In order to acquire ground truth measurements,
we propose using three total stations, each of them being labelled from one to
three. Each total station produces measurements in its own coordinate frame
$\mathcal{F}_{1},\mathcal{F}_{2}$ and $\mathcal{F}_{3}$. We also propose to
attach three prisms to the robotic platform, one per total station. These
prisms can be automatically detected by total stations in a _tracking mode_.
In this mode, each total station continuously tracks a single prism and
reports a measurements at a fixed rate. Also, a unique electronic signature of
each prism allows them to be distinguished while they are close to each other.

After installing the total stations on the experiment site, a calibration
procedure is required to find the rigid transformation
$\prescript{1}{2}{\mathbf{T}}$ from $\mathcal{F}_{2}$ to $\mathcal{F}_{1}$ and
the rigid transformation $\prescript{1}{3}{\mathbf{T}}$ from $\mathcal{F}_{3}$
to $\mathcal{F}_{1}$. This step is necessary to allow the expression of all
total station measurements in a common coordinate frame, which is the one of
total station 1. Small reflective targets are positioned around the site in
such a manner that they are visible to each of the total stations. The
measurements are done using the _single measurement_ mode, which requires
manually pointing the total stations to the targets.

Our setup requires that the data coming from all total stations be
synchronized according to one central clock. The clock chosen is the one on
the embedded computer located on the robotic platform. This conveniently
allows time synchronization with the robot’s own sensor measurements, while
providing it with real-time pose updates. However, because of the long range
tracking ability of the total stations, the distance between them and the
robot can be high. To solve this issue, long-range radio communication systems
are necessary.

We aim for real-time data gathering and broadcasting instead of having files
stored on the disk and post-processed later on. Real-time tracking can allow
the use of geo-fences for safety or the replacement of the on-board
localization algorithm when testing other algorithms. Each total station is
connected to its own embedded computer with a wired USB connection. These
embedded computers run a process which initiates the total station
measurements and receives the raw data. The raw data of the total station
measurements are the vertical angle, horizontal angle, range to the prism,
measurement time and the total station status. These data are then marshalled
by each of these embedded computers, dubbed _clients_ , via a long-range radio
module over to the _master_ embedded computer, located in the mobile robotic
platform. The purpose of this master is to synchronize the data coming from
the clients with a single clock and to control data flow between all long-
range radios via a remote data collection protocol. For broader acceptance and
portability, we restrict ourselves to single-channel communication. Because of
this, the master must poll one client at a time to collect the data. Our
communication protocol is composed of two parts: a synchronization protocol
between the master and the clients, and the sending of the data.

Figure 2: The experimental setup composed of three total stations, a Clearpath
Warthog with prisms (golden, blue, sienna) and two GPS receivers (pink). The
total station coordinate frames are related by the rigid transforms
$\prescript{1}{2}{\mathbf{T}}$ and $\prescript{1}{3}{\mathbf{T}}$. We find
these transforms by placing reflective targets (bottom right) around the
experimental area and then measuring their position by all three total
stations. The pose of of the robot is expressed by the transform
$\prescript{1}{R}{\mathbf{T}}$ from the robot frame $\mathcal{F}_{R}$ to
$\mathcal{F}_{1}$.

###  III-B Time synchronization

The system clocks on the clients are not synchronized with the system clock on
the master. The synchronization protocol is used to find the skew between the
clocks of the clients and the master. To begin a synchronization, the master
sends a message to ask a client to switch into this mode. One client at a time
can be synchronized. The protocol is composed of cycles using four messages
which is divided into two steps. The number of cycles are set according to the
configuration. The general timeline of the communication protocol is shown in
the top of Figure 3, and details about the synchronization procedure are in
the bottom of the same figure.

Figure 3: Communication protocol between the clients and the master radio
module, with the detail of the time synchronization protocol.

First, two concurrent time stamps $t_{m,i}$ and $t_{i}$ are captured in the
master and in the client. The subscript $i$ is the cycle index, $m$ stands for
_middle_. The two first messages are used for this purpose. At time $t_{b,i}$
($b$ as _beginning_), the master sends a very short message (two bytes) to the
client and saves the value of $t_{b,i}$. Upon receipt of the message, the
client stores the current time $t_{i}$ and also sends a very short reply to
confirm the reception of the request. The master receives this confirmation at
time $t_{e,i}$ ($e$ as _end_), and records the value of $t_{e,i}$. The
timestamp $t_{m,i}$ can be estimated for every cycle as the average of
$t_{b,i}$ and $t_{e,i}$.

In the second step, we compute the skew between the master clock and the
client clock. By the third message of the cycle, the master requests the value
$t_{i}$ from the client. The last message of the cycle is the reply of $t_{i}$
value back to the master. Using the value $t_{m,i}$, the skew of each cycle is
found by subtracting the times $t_{i}$ and $t_{m,i}$. Then, an average of the
skew $\overline{d}$ is calculated over all $i$ cycles. This average value is
the correction which is applied on the timestamp of the data incoming from the
corresponding client.

An initial synchronization is conducted at the beginning of the operation. All
of the clients are synchronized one at a time. Then, the estimated skews are
directly applied to correct the timestamps of the data. Since the clients have
a clock with significant drift, we need to synchronize them periodically
during the operation. To avoid losing too much measurement data, we only re-
synchronize one client at a time during an experiment. These re-
synchronizations are done following the same protocol as during the initial
synchronization, yet with less cycles. Because less messages are sent to
compute the time difference between the master and the client, the corrections
calculated are not directly applied to the data timestamp. A low-pass Infinite
Impulse Response (IIR) filter, which takes the former corrections
$\delta_{j-1}$ calculated into account, is applied using

$\delta_{j}=w\overline{d}+(1-w)\delta_{j-1}.$ ((1))

The final correction $\delta_{j}$ is applied after the $j^{th}$
synchronization. The filter weight $w$ is typically set $w\ll 0.5$ to slowly
update the time offset. This filter is used only after the initial
synchronization.

After the synchronization protocol, the master switches the client to the
normal measurement mode. Following this protocol, the total stations can
transmit the data, one total station at a time due to the single-channel
available. The master sends a measurement-request message to the total station
in turn. If a new measurement is available, it is sent back by the client.
Otherwise, the master passes to the next total station and this communication
cycle repeats itself.

###  III-C Calibration procedure

To find the rigid transformations $\prescript{1}{2}{\mathbf{T}}$ and
$\prescript{1}{3}{\mathbf{T}}$, a point-to-point alignment minimization is
carried out by using the Singular Value Decomposition (SVD) algorithm [24]
between the manually-measured reflective marker positions in the common frame
$\mathcal{F}_{1}$ and the other two total station frames, $\mathcal{F}_{2}$
and $\mathcal{F}_{3}$. Cartesian coordinates of each total station measurement
can be projected into $\mathcal{F}_{1}$ using $\prescript{1}{2}{\mathbf{T}}$
and $\prescript{1}{3}{\mathbf{T}}$.

It is also necessary to define the positions of the three prisms on the robot,
with respect to the robot coordinate frame. This measurement is performed
accurately by one total station only once when mounting the prisms on the
robot. With the help of a Computer-aided design (CAD) software and the model
of the robot, the prisms can be localized with respect to any chosen
coordinate frame.

###  III-D Six-DOF pose estimation from prism reflector position measurements

After calibration and synchronization of the system, we proceed to the actual
measurements. First, the validity of the data is checked. If a total station
reports a status indicating an error during the measurement (e.g., prism not
detected, prism too close, total station not levelled), the measurement is not
taken into account. Once checked, the raw measurements are transformed into
Cartesian coordinates in the corresponding total station frame and then into
the common $\mathcal{F}_{1}$ frame.

Once the prism position measurements all lie in the same frame, interpolation
is performed as shown in Figure 4. The figure shows the process for a single
axis $X$, but the same is done for $Y$ and $Z$. For a given time $t_{i}$, we
interpolate the position of each of the three prisms between the actual
measurements. In the case of measurement outage of any total station for a
time longer than a given threshold, the corresponding interval is marked as
invalid for obtaining the ground truth reference. An example of this situation
is depicted in the grey interval of Figure 4. The chosen time step for
interpolation is $50\text{\,}\mathrm{ms}$ to obtain a $20\text{\,}\mathrm{Hz}$
sampling rate in the reference. Results are triplets of coordinates
corresponding to the positions of the three prisms at time $t_{i}$.

Figure 4: Example of measurement data interpolation for a single axis. Each
color represents a single tracked prism. Solid shapes are actual measurements,
while hollow shapes are interpolated estimations at time $t_{i}$.

For each triplet, a point-to-point alignment minimization is performed, using
the same method as the one used during the
$\mathcal{F}_{1}$-$\mathcal{F}_{2}$-$\mathcal{F}_{3}$ calibration procedure.
The problem can be summarized as minimizing the point-to-point cost function
$\mathrm{J}_{\text{p-p}}(\cdot)$ as

$\displaystyle\prescript{1}{R}{\mathbf{T}}=\arg\min_{\bm{T}}\>{J}_{\text{p-p}}\left(\mathcal{Q},\mathcal{P}\right)\text{,
\hskip 1.42271pt with \hskip
1.42271pt}{J}_{\text{p-p}}=\sum_{k=1}(\bm{e}_{k}^{T}\bm{e}_{k})\text{,}$ ((2))

where $\bm{e}_{k}=\bm{q}_{k}-\bm{T}\bm{p}_{k}$ is the error vector between the
$k$th point $\bm{q}$ of $\mathcal{Q}$ representing the measured prism
positions taken as reference, and the $k$th point $\bm{p}$ of $\mathcal{P}$
representing the positions of the prisms as mounted on the robot and
calibrated in the lab, and $\bm{T}\in\mathrm{SE}(3)$ is a rigid transformation
matrix used during the minimization. The result given by (2) is a rigid
transformation matrix $\prescript{1}{R}{\mathbf{T}}\in\mathrm{SE}(3)$ between
the frame $\mathcal{F}_{1}$ and the frame of the robot $\mathcal{F}_{R}$,
which gives the 6-DOF pose of the robot body frame (a translation vector and a
rotation quaternion) in the common frame $\mathcal{F}_{1}$.

##  IV Experiments

To evaluate the performance of our position referencing system, we conducted
several experiments in different locations around the Montmorency Forest.
Later is a research facility belonging of Laval University, located north of
Quebec City. The specific locations were _i)_ one experiment in a former
quarry , _ii)_ three experiments in a forest cross-country ski trail, and
_iii)_ one experiment in a forest snowmobile trail. Figure 1 and Figure 5 show
these environments. These locations were chosen to test our setup in different
operating conditions: medium path of around $200\text{\,}\mathrm{m}$ in an
open space for the quarry, short path of around $150\text{\,}\mathrm{m}$ with
foliage and trees on the ski trail, and longer path of around
$600\text{\,}\mathrm{m}$ in the forest on the snowmobile trail. Hence, we
expect some GNSS signal attenuation and some long-range radio signal
attenuation in the forest area. These will allow us to evaluate our total
stations setup compared to the accuracy of a GNSS setup, but also to evaluate
the tracking mode of the total station over different experiments. The
evaluations will compare the results obtained during each of these
experiments.

Figure 5: Experimental environments. _Left_ : the ski trail. _Right_ : the
stone quarry areas.

###  IV-A Equipment used

The particular equipment used for our experiments is displayed in Figure 2. In
this figure, we can see a _ClearPath Warthog_ mobile robot used as the test
robotic platform, along with three _Trimble MultiTrack_ prisms on the
platform. These active prisms emit an infrared signature producing a unique
identifier to avoid confusion during tracking. They were individually tracked
by three total stations _Trimble S7_. In nominal conditions, the distance
measurement has an accuracy of $2\text{\,}\mathrm{mm}$ and the range of the
tracking mode is $800\text{\,}\mathrm{m}$. The angular accuracy is
$1\text{\,}\mathrm{{}^{\prime}}\mathrm{{}^{\prime}}$ for the vertical and
horizontal angles. This means that at $800\text{\,}\mathrm{m}$, the angular
accuracy will affect the precision by an error of almost
$4\text{\,}\mathrm{mm}$. The maximum achievable measurement rate of a _Trimble
S7_ in the prism-tracking mode is $2.5\text{\,}\mathrm{Hz}$.

The radio communication between all of the total stations and the master
computer is provided by four _Raspberry Pi 4_ computers with _LoRa Shield_
radio modules from _Dragino_. The parameters used for the communication are
listed in Table I. The chosen modulation allows the Raspberry Pi to reliably
send data over a distance of up to $800\text{\,}\mathrm{m}$ at 366 Bps in open
space. With our communication protocol, we can poll data from all three total
stations at a frequency of approximately $1.4\text{\,}\mathrm{Hz}$.

TABLE I: Table of the radio protocol communication parameters used. _Parameters_ | _Values_
---|---
LoRa frequency | $905\text{\,}\mathrm{MHz}$
Spreading factor | 9
Error coding rate | 4:6
Maximum byte rate | $366\text{\,}\mathrm{B}\mathrm{p}\mathrm{s}$
Measurement message length | $34\text{\,}\mathrm{B}$
Low-pass filter parameter $w$ | 0.1
Initial synchronization cycles | 50
Re-synchronization cycles | 5
Measurement outage threshold | $1\text{\,}\mathrm{s}$

Three _Reach RS+_ RTK receivers were used for comparison purposes, two on the
robot and one as a reference on a tripod. The distance between each prism, as
well as the distance between the two GNSS receivers on the robot were
accurately measured by a single total station in the lab. The obtained values
are shown in the Table II. To assess the accuracy of our localization
approach, we compare the reference values from Table II and the ones measured
during the experiments for noise evaluation, similarly to the work of [23].

TABLE II: Distance between the prisms and between the GNSS receivers, used as reference measurements. | _Prisms 1-2_ | _Prisms 1-3_ | _Prisms 2-3_ | _GPS 1-2_
---|---|---|---|---
_Distance_ |  $987\text{\,}\mathrm{mm}$ $\pm$ $2\text{\,}\mathrm{mm}$ |  $681\text{\,}\mathrm{mm}$ $\pm$ $2\text{\,}\mathrm{mm}$ |  $815\text{\,}\mathrm{mm}$ $\pm$ $2\text{\,}\mathrm{mm}$ |  $810\text{\,}\mathrm{mm}$ $\pm$ $2\text{\,}\mathrm{mm}$

###  IV-B Evaluation in static conditions

We first evaluate the precision of our system with the robot during stationary
poses from data taken at the ski trail shown in Figure 5-_left_. We identified
the intervals when the robot was stationary based on the recorded odometry.
For these intervals, the distances between the measured prism positions were
calculated. We used the raw data of the total stations without the
interpolation for this purpose. This was done to evaluate the accuracy of the
calibration procedure, which does not depend on the interpolation, and also to
evaluate the precision of the total stations themselves. The 345 resulting
distance samples are aggregated in the histograms shown in Figure 6. We can
compare the mean values $\mu_{12}$, $\mu_{13}$, and $\mu_{23}$ with the
reference values from the Table II, shown as the golden lines in the the
histograms. The absolute errors are respectively
$0.46\text{\,}\mathrm{m}\mathrm{m}$, $0.2\text{\,}\mathrm{m}\mathrm{m}$, and
$0.13\text{\,}\mathrm{m}\mathrm{m}$. The precision is high with the maximum
standard deviation of $3.5\text{\,}\mathrm{mm}$ in the case of $\sigma_{13}$.
This is in line with the expected precision of the total stations using the
tracking mode with a prism, given the $2\text{\,}\mathrm{mm}$ RMSE
specification from the manufacturer (_Trimble_). This also confirms that the
calibration procedure used to compute the common frame $\mathcal{F}_{1}$ was
accurate enough. We can also observe a symmetrical result around the true
value of the inter-prism distance 1 to 3, which results in a larger dispersion
$\sigma_{13}$ of the data than the $\sigma_{12}$ and $\sigma_{23}$ dispersion
of the data. We believe this is due to the calibration precision between the
total station frames, which caused a small error in rotation in the transform
matrix $\prescript{1}{3}{\mathbf{T}}$. This is due to the experiment itself,
where the robot moved towards and then away from the total stations.

Figure 6: Histogram of the distance between prisms during stationary poses.
The golden lines represent the reference distance between the prisms given in
Table II.

###  IV-C Evaluation in dynamic conditions

Next, we evaluate the performance during the motion of the robot. Figure 7
shows the inter-prism error over time obtained from one ski trail experiment.
During this experiment, the robot was driven in a start-and-stop motion
fashion. The plot shows a representative portion of the run, between
$80\text{\,}\mathrm{s}$ and $140\text{\,}\mathrm{s}$, as the rest of the run
was similar. The linear speed and angular speed of the robot during this
experiment are shown in the bottom part of the figure. These graphs show that
the reference is less precise when the robot starts to move, or when it stops.
At these moments, the error increases, leading to the several spikes in the
distance error plot. Also, there is an increase in the error when the robot
turns at time $109\text{\,}\mathrm{s}$. We attribute this error to the ability
of the total stations to accurately track the prisms during abrupt changes of
velocity. Also, the sampling scheme dictated by our hardware has an influence
on the precision in this dynamic case. One needs to recall that the system has
a relatively low sampling rate and the measurements of the three prisms are
not taken synchronously.

Figure 7: Results in dynamic conditions in the ski trail experiment. _Top_ :
the inter-prism errors. _Bottom_ : the corresponding linear and angular speed
of the robot.

To examine the influence of the linear speed, acceleration and angular speed,
all data taken during the three ski-trail experiments were gathered together.
To each inter-prism distance error, the values of linear speed, acceleration
and angular speed were attributed, as depicted in Figure 8. All of the samples
belonging to a single bin were averaged. The first subplot (a) represents the
inter-prism distance error according to the linear speed and the angular speed
of the robot. We see a higher error at linear speeds between
$0.25\text{\,}\mathrm{m}\mathrm{/}\mathrm{s}$ and
$0.5\text{\,}\mathrm{m}\mathrm{/}\mathrm{s}$. The error also increases with
the angular speed. We think that these results have a correlation with the way
a human operator drives the robot. The robot driven by the operator either
goes straight forward at higher velocity or turns in lower velocities. This
human bias is also why the angular velocities better represented in the plot
at low linear speeds. When moving straight forward, our sampling rate is
sufficient for the interpolation to return correct values. However, when the
robot turns, we can miss some key samples and thus obtain less precise
results. The second subplot (b) shows the inter-prism distance error as a
function of acceleration and angular speed. The error spikes due to the start-
stop motion shown in the Figure 7 are clearly highlighted at both extremes of
the acceleration axis.

(a)

(b)

Figure 8: The distance error between prisms during all ski trail experiments.

###  IV-D From prism uncertainties to pose uncertainty

We analyze the impact of uncertainty in the prism measurements on the measured
robot position and orientation. As explained in Section III-D, after the
interpolation, a point-to-point alignment minimization is performed to obtain
the 6-DOF pose of the robot frame. To determine the level of uncertainty in
the 6-DOF pose that results from the noise in the total station measurements,
we carried out the following perturbation analysis: a Gaussian white noise was
randomly applied to the positions of the prisms with the standard deviation
range of $0\text{\,}\mathrm{mm}$ to $400\text{\,}\mathrm{mm}$ with a step of
$4\text{\,}\mathrm{mm}$. This approach was used to simplify the study of the
noise propagation through the point-to-point minimization, which uses non-
linear functions. The range used is comparable to the uncertainty that is
observed during our experiments. For each value in the range of standard
deviation, 1000 simulations of point-to-point alignment minimization were
performed. The resulting sets of 1000 poses allowed us to study the
distribution of position and orientation errors. Figure 9 shows these pose
errors as a function of the input perturbation. The mean and standard
deviation are extracted and plotted for each axis in position and orientation.
In this analysis, we choose the Euler angles representation of the orientation
in Figure 9 for the simplicity of their interpretation. The values are fairly
low and thus far from the singularities in this representation. Not
surprisingly, the spread of the inter-prism error position is very close to a
linear function depending on the noise uncertainty. The same applies for the
Euler angles, at least in the noise range we have chosen. Interestingly, the
mean value is not zero. Its absolute value is slightly increasing as a non-
linear function of noise range. Nevertheless in our experiments, the mean
value of the inter-prism measurement error is $10\text{\,}\mathrm{mm}$ and the
vast majority of its distribution lies below $50\text{\,}\mathrm{mm}$. In this
region, the observed bias in the resulting position and orientation can be
neglected. Moreover, the mean value of the inter-prism measurement error of
$10\text{\,}\mathrm{mm}$ translates to $10\text{\,}\mathrm{mm}$ in position
and $0.01\text{\,}\mathrm{r}\mathrm{a}\mathrm{d}$ ($0.6\text{\,}\deg$) in
orientation, based on sampling from an equivalent input error distribution.

Figure 9: Perturbation analysis of the 6-DOF pose estimation from the prism
measurements. The solid line and the shaded area represent respectively the
mean value and the standard deviation.

###  IV-E Comparison with GNSS solutions

To ground the performance of the total station referencing system to a
commonly-known alternative, we compared it to the RTK solution. In our
previous work [8], we have studied the impact of vegetation and obstacles both
on the propagation of the satellite signal and the base station correction
signal. The satellite localization performed well in open areas, but a
decrease in accuracy was observed when portions of the sky view were
obstructed by vegetation. To this effect, we have recorded the RTK positioning
data during two experiments: one in the quarry and another one in the forest.
The raw GNSS data were recorded as log files in the receivers and the RTK
solution was obtained in post-processing by the means of the _RTKLIB_
111http://www.rtklib.com/ library. We selected the RTK solver settings that
would be achievable in real-time operation (forward-pass filtering only).

Similarly to using the inter-prism distance to estimate the total station
localization precision, we rely on the measured distance between the GNSS
receivers to estimate their precision. The distance between these GNSS
receivers is indicated in Table II. Figure 10 shows the results from the
quarry experiment and from the snowmobile trail experiment. The purpose of
these experiments is to compare the results coming from an open-area and from
a more obstructed one. The average number of satellites used by the receivers
during the quarry experiment was $15.3\pm 1.4$ (minimum 2 and maximum 20
satellites), and $13.1\pm 2.8$ satellites in the forest area during the
snowmobile trail experiment (minimum 3 and maximum 21 satellites). For
reference, we show the precision of the basic GNSS solution without the RTK
corrections in Figure 10 (a). Even in the open area of a quarry, the precision
is essentially unable to accurately locate the robot. In a more realistic
comparison, Figure 10 (b) compares the precision of the corrected RTK
localization and our total station referencing system. In this environment and
with the robot moving, the precision is equivalent as shown on the histogram.
The total station system has a slightly higher error spread compared to the
RTK solution due to the error spike in the beginning of the experiment. This
error arises from an incorrect inter-prism distance 1-2, probably due to the
acceleration of the robot when it begins to move, similarly as in Figure 7.
Otherwise, the accuracy is similar for both. In a more difficult scenario,
Figure 10 (c) shows portions of the snowmobile trail experiment, where the
robot was in reach of all three total stations. From the beginning until the
mark of $200\text{\,}\mathrm{s}$ and in the remainder of the experiment (right
plot), the GNSS signals are significantly attenuated by the trees growing on
both sides of the trail. During this part, the GNSS precision drops
drastically, hovering around $1.6\text{\,}\mathrm{m}$. From the mark
$200\text{\,}\mathrm{s}$, the sky view opens and the GNSS achieves its nominal
precision again. On the other hand, the average precision of the total station
reference is $7.0\text{\,}\mathrm{mm}$ as long as the robot lies in the field
of view of the total stations.

(a)

(b)

(c)

Figure 10: The distance errors of our system compared to GPS during the quarry
and the forest experiments.

These two experiments show that the proposed total station setup has an
equivalent accuracy and precision compared to the RTK solution in an open
area. However, our solution is significantly better in an obstructed
environment such as a forest, where the GNSS signal is weaker because of the
vegetation. Our setup can thus be used as a source of precise ground truth in
such an environment.

##  V Conclusion

In this paper, we proposed a referencing system based on total stations in
order to provide 6-DOF pose measurements in real time. The referencing system
was demonstrated successfully on a large tracked skid-steered robot in a
subarctic forest. With the robot not moving, the precision achieved was better
than $1\text{\,}\mathrm{mm}$ in position and $0.06\text{\,}\deg$ in
orientation. These values are in line with the limits of the total stations
used for the experiments. With the robot moving, the precision decreased to
$10\text{\,}\mathrm{mm}$ in position and $0.6\text{\,}\deg$ in orientation. In
this dynamic case, the other factors are the ability of the total stations to
quickly react to speed changes and the measurement rate, which may be hiding
the dynamics. Nevertheless, the experiments have demonstrated the robustness
of our results against different types of environments. Also, they are
equivalent to the results obtained by the RTK solution in open areas and by a
large margin better in environments where GNSS signals are weak.

In future work, the reflective calibration markers will be replaced by a
dedicated, automatically tracked prism. This will increase the inter-total
station calibration accuracy and significantly decrease the setup time. We
will study in more detail the effect of acceleration of the robot, which we
consider to be the main source of error in our experiments at the moment. A
similar problem is addressed by [25], where they propose improvements based on
data from an RGB camera. This could be a starting point to address this
limitation. The total stations can achieve the maximum measurement rate of
$2.5\text{\,}\mathrm{H}\mathrm{z}$. We will further optimize the long-range
communication protocol to fully use this (limited) potential. Finally, we have
to assess the influence of prism configuration on the robot on the uncertainty
of the final pose, and possibly finding an optimal configuration.

## Acknowledgment

This research was supported by the Natural Sciences and Engineering Research
Council of Canada (NSERC) through the grant CRDPJ 527642-18 SNOW (Self-driving
Navigation Optimized for Winter).

## References

  * [1] J. Delmerico and D. Scaramuzza  “A Benchmark Comparison of Monocular Visual-Inertial Odometry Algorithms for Flying Robots”  In _2018 IEEE International Conference on Robotics and Automation (ICRA)_ , 2018, pp. 2502–2509  DOI: 10.1109/ICRA.2018.8460664
  * [2] Pierre Merriaux et al.  “A Study of Vicon System Positioning Performance”  In _Sensors_ 17.7, 2017, pp. 1–18  DOI: 10.3390/s17071591
  * [3] Joshua S. Furtado et al.  “Comparative Analysis of OptiTrack Motion Capture Systems”  In _Advances in Motion Sensing and Control for Robotic Applications_ , 2019, pp. 15–31
  * [4] Lingfeng Sang, Hongbo Wang, Hongnian Yu and Luige Vladareanu  “A novel human-carrying quadruped walking robot”  In _International Journal of Advanced Robotic Systems_ 14.4, 2017, pp. 1–17  DOI: 10.1177/1729881417716592
  * [5] Y. Morales and T. Tsubouchi  “DGPS, RTK-GPS and StarFire DGPS Performance Under Tree Shading Environments”  In _IEEE International Conference on Integration Technology_ , 2007, pp. 519–524  DOI: 10.1109/ICITECHNOLOGY.2007.4290370
  * [6] P McGarey et al.  “Field Report: Developing and Deploying a Tethered Robot to Map Extremely Steep Terrain”  In _Journal of Field Robotics_ 35.8, 2018, pp. 1241–1327
  * [7] Zichao Zhang, H. Rebecq, C. Forster and D. Scaramuzza  “Benefit of large field-of-view cameras for visual odometry”  In _IEEE International Conference on Robotics and Automation (ICRA)_ , 2016, pp. 801–808  DOI: 10.1109/ICRA.2016.7487210
  * [8] Vladimír Kubelka et al.  “Radio propagation models for differential GNSS based on dense point clouds”  In _Journal of Field Robotics_ 37.8, 2020, pp. 1347–1362  DOI: https://doi.org/10.1002/rob.21988
  * [9] P Krüsi et al.  “Lighting-invariant adaptive route following using iterative closest point matching”  In _Journal of Field Robotics_ 32.4, 2015, pp. 534–564
  * [10] Kirk MacTavish, Michael Paton and Timothy D. Barfoot  “Selective memory: Recalling relevant experience for long-term visual localization”  In _Journal of Field Robotics_ 35.8, 2018, pp. 1265–1292  DOI: 10.1002/rob.21838
  * [11] Patrick McGarey, Kirk MacTavish, François Pomerleau and Timothy D Barfoot  “TSLAM: Tethered simultaneous localization and mapping for mobile robots”  In _The International Journal of Robotics Research_ 36.12, 2017, pp. 1363–1386  DOI: 10.1177/0278364917732639
  * [12] David Wettergreen et al.  “Design and field experimentation of a prototype Lunar prospector”  In _International Journal of Robotics Research_ 29.12, 2010, pp. 1550–1564  DOI: 10.1177/0278364910370217
  * [13] Daniel M. Helmick et al.  “Path following Using visual odometry for a Mars rover in high-slip environments”  In _IEEE Aerospace Conference Proceedings_ 2, 2004, pp. 773–789  DOI: 10.1109/AERO.2004.1367679
  * [14] G Hitz, F Pomerleau, F Colas and R Siegwart  “Relaxing the planar assumption: 3D state estimation for an autonomous surface vessel”  In _International Journal of Robotics Research_ 34.13, 2015, pp. 1604–1621
  * [15] Marko Bjelonic et al.  “Weaver: Hexapod robot for autonomous navigation on unstructured terrain”  In _Journal of Field Robotics_ 35.7, 2018, pp. 1063–1079  DOI: 10.1002/rob.21795
  * [16] Vladimír Kubelka et al.  “Robust Data Fusion of Multimodal Sensory Information for Mobile Robots”  In _Journal of Field Robotics_ 32.4, 2015, pp. 447–473  DOI: 10.1002/rob.21535
  * [17] Daniel Loret Mola Lemus, David Kohanbash, Scott Moreland and David Wettergreen  “Slope Descent using Plowing to Minimize Slip for Planetary Rovers”  In _Journal of Field Robotics_ 31.5, 2014, pp. 803–819  DOI: 10.1002/rob.21518
  * [18] Juraj Kabzan et al.  “AMZ Driverless: The full autonomous racing system”  In _Journal of Field Robotics_ 37.7, 2020, pp. 1267–1294  DOI: 10.1002/rob.21977
  * [19] Michael Burri et al.  “The EuRoC micro aerial vehicle datasets”  In _International Journal of Robotics Research_ 35.10, 2016, pp. 1157–1163  DOI: 10.1177/0278364915620033
  * [20] Patrik Schmuck and Margarita Chli  “CCM-SLAM: Robust and efficient centralized collaborative monocular simultaneous localization and mapping for robotic teams”  In _Journal of Field Robotics_ 36.4, 2019, pp. 763–781  DOI: 10.1002/rob.21854
  * [21] Giuliano Punzo et al.  “Bipartite Guidance, Navigation and Control Architecture for Autonomous Aerial Inspections Under Safety Constraints”  In _Journal of Intelligent and Robotic Systems: Theory and Applications_ 95.3-4 Journal of Intelligent & Robotic Systems, 2019, pp. 1049–1061  DOI: 10.1007/s10846-018-0780-1
  * [22] T.. White et al.  “A mobile climbing robot for high precision manufacture and inspection of aerostructures”  In _International Journal of Robotics Research_ 24.7, 2005, pp. 589–598  DOI: 10.1177/0278364905055701
  * [23] François Pomerleau, Ming Liu, Francis Colas and Roland Siegwart  “Challenging data sets for point cloud registration algorithms”  In _The International Journal of Robotics Research_ 31.14, 2012, pp. 1705–1711  DOI: 10.1177/0278364912458814
  * [24] K.. Arun, T.. Huang and S.. Blostein  “Least-Squares fitting of two 3-D Point Sets”  In _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 1987, pp. 698–700
  * [25] Matthias Ehrhart and Werner Lienhart  “Object tracking with robotic total stations: Current technologies and improvements based on image data”  In _Journal of Applied Geodesy_ 11.3, 2017, pp. 131–142  DOI: 10.1515/jag-2016-0043

  *[DOFs]: Degrees Of Freedom
  *[RTK]: Real Time Kinematics
  *[GNSS]: Global Navigation Satellite System
  *[DOF]: Degrees Of Freedom
  *[IIR]: Infinite Impulse Response
  *[SVD]: Singular Value Decomposition
  *[CAD]: Computer-aided design


# Indirect Measurement of Switch Terms of a Vector Network Analyzer with
Reciprocal Devices

Ziad Hatab, Michael Ernst Gadringer, and Wolfgang Bösch
Graz University of Technology, Austria
{z.hatab; michael.gadringer<EMAIL_ADDRESS>Software code and measurements
are available online:
https://github.com/ZiadHatab/vna-switch-terms

###### Abstract

This paper presents an indirect method for measuring the switch terms of a
vector network analyzer (VNA) using at least three reciprocal devices, which
do not need to be characterized beforehand. This method is particularly
suitable for VNAs that use a three-sampler architecture, which allows for
applying first-tier calibration methods based on the error box model. The
proposed method was experimentally verified by comparing directly and
indirectly measured switch terms and performing a multiline thru-reflect-line
(TRL) calibration.

###### Index Terms:

VNA, calibration, microwave measurement

##  I Introduction

Calibration of a vector network analyzer (VNA) is crucial for removing
systematic errors between the measured device under test (DUT) and the actual
receivers of the VNA. The most common calibration method is the short-open-
load-thru (SOLT) method, which is based on the 12-term error model of a two-
port VNA. However, this method requires all standards to be fully
characterized. In [1], a modification to the SOLT method was introduced where
the thru standard is replaced with any transmissive reciprocal device, called
the SOLR method. This method is based on the error box model of a two-port
VNA. Other advanced self-calibration methods, including thru-reflect-line
(TRL), multiline TRL, line-reflect-match (LRM), and line-reflect-reflect-match
(LRRM) [2, 3, 4, 5], also rely on the error box model of a two-port VNA.

A limitation of calibration methods based on the error box model is that it
requires a four-sampler VNA to sample all waves, whereas the 12-term error
model can still be used in three-sampler VNAs. Fig. 1 illustrates the two VNA
sampling architectures. The difference between the two sampling architectures
is that in the three-sampler VNA, we do not sample the reflected wave of the
termination load of the non-driving port. Although the termination load is
generally designed to be matched, in reality, there is always some reflection
that needs to be accounted for. This reflection is called the switch term.
Since the ports are driven in both forward and reverse directions, there are
two switch terms.

(a)

(b)

Figure 1: Illustration of three- (a) and four-sampler (b) architectures of a
VNA. Both diagrams depict driving in the forward direction.

Since the terminations of the non-driving ports remain constant, and the
switching between the driving ports is often very repeatable with the help of
electronic switches, the switch terms introduce a systematic deviation and
only need to be measured once to be considered appropriately. These terms can
be regarded as part of the calibration coefficients by the conversion
relationships between the 12-term and error box models [6, 7, 8, 9].

In three-sampler VNAs, self-calibration methods based on the error box model
cannot be used as a first-tier calibration, as we cannot directly measure the
switch terms. However, error box calibration methods can be performed as a
second-tier calibration after a SOLT calibration [10]. Such methods require
pre-characterized calibration standards, which goes against the purpose of
self-calibration methods that use partially defined standards.

This paper aims to introduce a new method to indirectly measure the switch
terms using at least three transmissive reciprocal devices, which do not need
to be characterized beforehand. The proposed method enables the usage of error
box calibration methods in three-sampler VNAs without requiring any first-tier
calibration.

##  II Mathematical Formulation

###  II-A Problem statement

In a two-port VNA, when all four waves are sampled in both driving directions,
the measured S-parameters are described using the following notation [11]:

$\begin{bmatrix}\hat{b}_{11}&\hat{b}_{12}\\\
\hat{b}_{21}&\hat{b}_{22}\end{bmatrix}=\boldsymbol{S}\begin{bmatrix}\hat{a}_{11}&\hat{a}_{12}\\\
\hat{a}_{21}&\hat{a}_{22}\end{bmatrix}$ (1)

where $\hat{a}_{ij}$ and $\hat{b}_{ij}$ represent the sampled incident and
reflected waves, respectively, at port-$i$ when driven by port-$j$.

In a three-sampler VNA, the waves $\hat{a}_{12}$ and $\hat{a}_{21}$ are not
measured due to a lack of dedicated receivers. To address this, the measured
incident waves in (1) can be split into two matrices as follows:

$\begin{bmatrix}\hat{b}_{11}&\hat{b}_{12}\\\
\hat{b}_{21}&\hat{b}_{22}\end{bmatrix}=\boldsymbol{S}\begin{bmatrix}1&\frac{\hat{a}_{12}}{\hat{a}_{22}}\\\
\frac{\hat{a}_{21}}{\hat{a}_{11}}&1\end{bmatrix}\begin{bmatrix}\hat{a}_{11}&0\\\
0&\hat{a}_{22}\end{bmatrix}$ (2)

By taking the inverse of the diagonal matrix on the right-hand side of (2), we
obtain the conventionally measured ratios.

$\begin{bmatrix}\frac{\hat{b}_{11}}{\hat{a}_{11}}&\frac{\hat{b}_{12}}{\hat{a}_{22}}\\\
\frac{\hat{b}_{21}}{\hat{a}_{11}}&\frac{\hat{b}_{22}}{\hat{a}_{22}}\end{bmatrix}=\boldsymbol{S}\begin{bmatrix}1&\frac{\hat{a}_{12}}{\hat{a}_{22}}\\\
\frac{\hat{a}_{21}}{\hat{a}_{11}}&1\end{bmatrix}$ (3)

If we define the ratios on the left-hand side of (3) as the measured
S-parameters, we can then rewrite the remaining ratios on the right-hand side
as follows:

$\begin{bmatrix}\mkern 1.5mu\overline{\mkern-1.5muS\mkern-1.5mu}\mkern
1.5mu_{11}&\mkern 1.5mu\overline{\mkern-1.5muS\mkern-1.5mu}\mkern
1.5mu_{12}\\\ \mkern 1.5mu\overline{\mkern-1.5muS\mkern-1.5mu}\mkern
1.5mu_{21}&\mkern 1.5mu\overline{\mkern-1.5muS\mkern-1.5mu}\mkern
1.5mu_{22}\end{bmatrix}=\boldsymbol{S}\begin{bmatrix}1&\mkern
1.5mu\overline{\mkern-1.5muS\mkern-1.5mu}\mkern 1.5mu_{12}\Gamma_{12}\\\
\mkern 1.5mu\overline{\mkern-1.5muS\mkern-1.5mu}\mkern
1.5mu_{21}\Gamma_{21}&1\end{bmatrix}$ (4)

where $\overline{S}_{ij}$ represents the measured S-parameters and
$\Gamma_{ij}$ represents the switch terms of the VNA:

$\mkern 1.5mu\overline{\mkern-1.5muS\mkern-1.5mu}\mkern
1.5mu_{ij}=\frac{\hat{b}_{ij}}{\hat{a}_{jj}},\qquad\Gamma_{ij}=\frac{\hat{a}_{ij}}{\hat{b}_{ij}}$
(5)

The switch terms are formed by the ratios of the receivers of the non-driving
port. Therefore, they are independent of the measured DUT, as any influence
introduced by the DUT will be seen equally by both waves $\hat{a}_{ij}$ and
$\hat{b}_{ij}$. In general, the switch term corrected S-parameters are given
as follows:

$\boldsymbol{S}=\begin{bmatrix}\mkern
1.5mu\overline{\mkern-1.5muS\mkern-1.5mu}\mkern 1.5mu_{11}&\mkern
1.5mu\overline{\mkern-1.5muS\mkern-1.5mu}\mkern 1.5mu_{12}\\\ \mkern
1.5mu\overline{\mkern-1.5muS\mkern-1.5mu}\mkern 1.5mu_{21}&\mkern
1.5mu\overline{\mkern-1.5muS\mkern-1.5mu}\mkern
1.5mu_{22}\end{bmatrix}\begin{bmatrix}1&\mkern
1.5mu\overline{\mkern-1.5muS\mkern-1.5mu}\mkern 1.5mu_{12}\Gamma_{12}\\\
\mkern 1.5mu\overline{\mkern-1.5muS\mkern-1.5mu}\mkern
1.5mu_{21}\Gamma_{21}&1\end{bmatrix}^{-1}$ (6)

In the special case where the measured two-port device is transmissionless,
the switch terms $\Gamma_{ij}$ do not influence the measurements as
$\overline{S}_{21}=\overline{S}_{12}=0$.

Using a four-sampler VNA, we can directly measure $\Gamma_{ij}$ by connecting
any transmissive device and calculating the ratios according to the definition
in (5). Furthermore, we can measure the S-parameters directly using (1). In
contrast to the four-sampler VNA, a three-sampler one can only measure
$\overline{S}_{ij}$. Therefore, it is advantageous for three-sampler VNAs to
find a way to measure $\Gamma_{ij}$ without measuring the waves $\hat{a}_{12}$
and $\hat{a}_{21}$, and using only $\overline{S}_{ij}$ measurements.

###  II-B Proposed indirect measurement of the switch terms

Fig. 2 shows the error box model of a two-port VNA. Using T-parameters, the
measured DUT is given in terms of wave-parameter as follows [11]:

$\begin{bmatrix}\hat{a}_{11}&\hat{a}_{12}\\\
\hat{b}_{11}&\hat{b}_{12}\end{bmatrix}=\boldsymbol{E}_{\mathrm{L}}\boldsymbol{T}_{\mathrm{D}}\boldsymbol{E}_{\mathrm{R}}\begin{bmatrix}\hat{a}_{21}&\hat{a}_{22}\\\
\hat{b}_{21}&\hat{b}_{22}\end{bmatrix}$ (7)

where $\boldsymbol{E}_{\mathrm{L}}$ and $\boldsymbol{E}_{\mathrm{R}}$ are the
left and right error boxes, and $\boldsymbol{T}_{\mathrm{D}}$ is the actual
DUT.

Figure 2: Two-port VNA error box model.

We split the wave-parameter matrices in (7) into two matrices as follows:

$\begin{bmatrix}1&\frac{\hat{a}_{12}}{\hat{b}_{12}}\\\
\frac{\hat{b}_{11}}{\hat{a}_{11}}&1\end{bmatrix}\begin{bmatrix}\hat{a}_{11}&0\\\
0&\hat{b}_{12}\end{bmatrix}=\boldsymbol{E}_{\mathrm{L}}\boldsymbol{T}_{\mathrm{D}}\boldsymbol{E}_{\mathrm{R}}\begin{bmatrix}\frac{\hat{a}_{21}}{\hat{b}_{21}}&1\\\
1&\frac{\hat{b}_{22}}{\hat{a}_{22}}\end{bmatrix}\begin{bmatrix}\hat{b}_{21}&0\\\
0&\hat{a}_{22}\end{bmatrix}$ (8)

The above expression can be simplified by multiplying the inverse of the
diagonal matrix at the right-hand side. This step reduces all wave parameters
into ratios as follows:

$\begin{bmatrix}1&\frac{\hat{a}_{12}}{\hat{b}_{12}}\\\
\frac{\hat{b}_{11}}{\hat{a}_{11}}&1\end{bmatrix}\begin{bmatrix}\frac{\hat{a}_{11}}{\hat{b}_{21}}&0\\\
0&\frac{\hat{b}_{12}}{\hat{a}_{22}}\end{bmatrix}=\boldsymbol{E}_{\mathrm{L}}\boldsymbol{T}_{\mathrm{D}}\boldsymbol{E}_{\mathrm{R}}\begin{bmatrix}\frac{\hat{a}_{21}}{\hat{b}_{21}}&1\\\
1&\frac{\hat{b}_{22}}{\hat{a}_{22}}\end{bmatrix}$ (9)

The final simplification is to replace the ratios with the definitions
established in (5). The rearranged expression is presented in (10).

$\begin{bmatrix}1&\Gamma_{12}\\\ \mkern
1.5mu\overline{\mkern-1.5muS\mkern-1.5mu}\mkern
1.5mu_{11}&1\end{bmatrix}\begin{bmatrix}1/\mkern
1.5mu\overline{\mkern-1.5muS\mkern-1.5mu}\mkern 1.5mu_{21}&0\\\ 0&\mkern
1.5mu\overline{\mkern-1.5muS\mkern-1.5mu}\mkern
1.5mu_{12}\end{bmatrix}=\boldsymbol{E}_{\mathrm{L}}\boldsymbol{T}_{\mathrm{D}}\boldsymbol{E}_{\mathrm{R}}\begin{bmatrix}\Gamma_{21}&1\\\
1&\mkern 1.5mu\overline{\mkern-1.5muS\mkern-1.5mu}\mkern
1.5mu_{22}\end{bmatrix}$ (10)

Our goal is to extract $\Gamma_{21}$ and $\Gamma_{12}$ without prior knowledge
of the error boxes or the DUT. We can do this by assuming that the DUT is a
reciprocal device, i.e.,
$\mathrm{det}\left(\boldsymbol{T}_{\mathrm{D}}\right)=1$. By applying the
determinate operator to (10) and using the property that
$\mathrm{det}\left(\boldsymbol{A}\boldsymbol{B}\right)=\mathrm{det}\left(\boldsymbol{A}\right)\mathrm{det}\left(\boldsymbol{B}\right)$,
we can derive the following:

$(1-\mkern 1.5mu\overline{\mkern-1.5muS\mkern-1.5mu}\mkern
1.5mu_{11}\Gamma_{12})\frac{\mkern
1.5mu\overline{\mkern-1.5muS\mkern-1.5mu}\mkern 1.5mu_{12}}{\mkern
1.5mu\overline{\mkern-1.5muS\mkern-1.5mu}\mkern
1.5mu_{21}}=\underbrace{\mathrm{det}\left(\boldsymbol{E}_{\mathrm{L}}\right)\mathrm{det}\left(\boldsymbol{E}_{\mathrm{R}}\right)}_{=c\text{
(constant)}}(\Gamma_{21}\mkern 1.5mu\overline{\mkern-1.5muS\mkern-1.5mu}\mkern
1.5mu_{22}-1)$ (11)

The above expression can be simplified as follows:

$\frac{\mkern 1.5mu\overline{\mkern-1.5muS\mkern-1.5mu}\mkern
1.5mu_{12}}{\mkern 1.5mu\overline{\mkern-1.5muS\mkern-1.5mu}\mkern
1.5mu_{21}}-\mkern 1.5mu\overline{\mkern-1.5muS\mkern-1.5mu}\mkern
1.5mu_{11}\frac{\mkern 1.5mu\overline{\mkern-1.5muS\mkern-1.5mu}\mkern
1.5mu_{12}}{\mkern 1.5mu\overline{\mkern-1.5muS\mkern-1.5mu}\mkern
1.5mu_{21}}\Gamma_{12}-\mkern 1.5mu\overline{\mkern-1.5muS\mkern-1.5mu}\mkern
1.5mu_{22}c\Gamma_{21}+c=0$ (12)

From (12), we can recognize that we have a linear equation in three unknowns:
$\Gamma_{12}$, $c\Gamma_{21}$, and $c$. Therefore, if we measure at least
three unique transmissive reciprocal devices, we can solve for these unknowns
by solving the following linear system of equations:

$\begin{bmatrix}-\mkern 1.5mu\overline{\mkern-1.5muS\mkern-1.5mu}\mkern
1.5mu_{11}^{(1)}\frac{\mkern 1.5mu\overline{\mkern-1.5muS\mkern-1.5mu}\mkern
1.5mu_{12}^{(1)}}{\mkern 1.5mu\overline{\mkern-1.5muS\mkern-1.5mu}\mkern
1.5mu_{21}^{(1)}}&-\mkern 1.5mu\overline{\mkern-1.5muS\mkern-1.5mu}\mkern
1.5mu_{22}^{(1)}&1&\frac{\mkern
1.5mu\overline{\mkern-1.5muS\mkern-1.5mu}\mkern 1.5mu_{12}^{(1)}}{\mkern
1.5mu\overline{\mkern-1.5muS\mkern-1.5mu}\mkern 1.5mu_{21}^{(1)}}\\\
\vdots&\vdots&\vdots&\vdots\\\ -\mkern
1.5mu\overline{\mkern-1.5muS\mkern-1.5mu}\mkern 1.5mu_{11}^{(M)}\frac{\mkern
1.5mu\overline{\mkern-1.5muS\mkern-1.5mu}\mkern 1.5mu_{12}^{(M)}}{\mkern
1.5mu\overline{\mkern-1.5muS\mkern-1.5mu}\mkern 1.5mu_{21}^{(M)}}&-\mkern
1.5mu\overline{\mkern-1.5muS\mkern-1.5mu}\mkern
1.5mu_{22}^{(M)}&1&\frac{\mkern
1.5mu\overline{\mkern-1.5muS\mkern-1.5mu}\mkern 1.5mu_{12}^{(M)}}{\mkern
1.5mu\overline{\mkern-1.5muS\mkern-1.5mu}\mkern
1.5mu_{21}^{(M)}}\end{bmatrix}\begin{bmatrix}\Gamma_{12}\\\ c\Gamma_{21}\\\
c\\\ 1\end{bmatrix}=\boldsymbol{0}$ (13)

where $M\geq 3$ is the number of measured reciprocal devices. We need at least
three distinct measurements to determine the unknowns, since the system matrix
must have a rank of 3 to be solved. In general, the uniqueness of the
reciprocal standards affects the conditioning of the system matrix.

To solve for the unknowns, we need to find the nullspace of the system matrix.
We can estimate the nullspace by applying the singular value decomposition
(SVD) to (13). In this case, the best approximate of the nullspace corresponds
to the right singular vector that is associated to the smallest singular value
[12]. Since the nullspace is only unique up to a scalar multiple, we can solve
for the switch terms by taking the ratio of the elements of the nullspace
vector as follows:

$\Gamma_{12}=\frac{v_{41}}{v_{44}},\qquad\Gamma_{21}=\frac{v_{42}}{v_{43}}$
(14)

where $\boldsymbol{v}_{4}=[v_{41},v_{42},v_{43},v_{44}]^{T}$ is the nullspace
vector found through the SVD.

##  III Experiment

The experiment consists of two parts. In the first step, we tested our
proposed method for extracting switch terms using only three of the four
available receivers on a VNA. We compared the obtained results to the switch
terms computed directly using the fourth receiver. As the second step, we
performed a multiline TRL calibration using both switch terms computed
directly and indirectly. The results were compared by calibrating a stepped
impedance line.

The reciprocal devices that we used consisted of a line standard from the
multiline TRL kit (50 mm line) and a series-shunt (L-circuit) of $100\,\Omega$
resistors, which were measured twice by flipping the ports, as it is an
asymmetric device ($S_{11}\neq S_{22}$). For the multiline TRL kit, we
implemented microstrip lines on an FR4 substrate with a trace width of 3 mm
and a substrate height of 1.55 mm. The lengths of the lines (referenced to the
first line) are as follows: $\\{0,2.5,10,15,50\\}\,\mathrm{mm}$, with the
reflect standard implemented as a short. The standards are shown in Fig. 3.

(a)

(b)

(c)

Figure 3: Measured structures. (a) microstrip line multiline TRL kit
($50\,\Omega$), (b) stepped impedance line ($90\,\Omega$), and (c) series-
shunt $100\,\Omega$ circuit.

The R&S ZVA is the four-sampler VNA used in this experiment. To extract the
switch terms, we measured the aforementioned reciprocal devices and processed
the data offline with the help of the scikit-rf package in Python [13]. The
results, along with those obtained by direct wave ratio computation with the
fourth receiver, are shown in Fig. 4. The presented results highlight that
both results overlap in magnitude and phase.

Figure 4: Comparison of direct and indirect measurements of the switch terms.

Finally, we performed a multiline TRL calibration using the algorithm in
reference [14]. In Fig. 5, we present the calibrated results of a stepped
impedance line in different scenarios: ignoring the switch terms
($\Gamma_{21}=\Gamma_{12}=0$), directly measuring the switch terms, and
indirectly measuring the switch terms. The results in Fig. 5 show that both
directly and indirectly measured switch terms deliver the same results.
However, ignoring the switch terms results in noise-like behavior on the
traces.

Figure 5: Results of the multiline TRL calibrated stepped impedance line.

##  IV Conclusion

In this paper, we presented an innovative method that utilizes a minimum of
three transmissive reciprocal devices to measure the switch terms of a VNA,
relying on just three receivers. We applied our method to compute the switch
terms using only three receivers through practical measurements conducted on a
four-sampler VNA. To validate the accuracy of our results, we compared them
with direct measurements obtained using the fourth receiver, revealing
comparable outcomes. Furthermore, we successfully demonstrated the
implementation of a first-tier multiline TRL calibration, utilizing only three
receivers of the VNA.

One significant advantage of the proposed method is that it does not require
prior knowledge of the reciprocal devices. For example, one could use
electronically controlled resistors at the test ports to quickly obtain the
switch terms. This approach is particularly useful for error box calibration
methods in multiport VNAs. By using this method, $N+1$ samplers can replace a
full-reflectometry architecture with $2N$ samplers, where $N$ is the number of
ports. Such a simplification in the setup of a multiport VNA significantly
reduces the complexity and cost of the device.

## Acknowledgment

The financial support by the Austrian Federal Ministry for Digital and
Economic Affairs and the National Foundation for Research, Technology, and
Development is gratefully acknowledged.

## References

  * [1] A. Ferrero and U. Pisani, “Two-port network analyzer calibration using an unknown ’thru’,” _IEEE Microwave and Guided Wave Letters_ , vol. 2, no. 12, pp. 505–507, 1992, doi: 10.1109/75.173410.
  * [2] G. Engen and C. Hoer, “Thru-reflect-line: An improved technique for calibrating the dual six-port automatic network analyzer,” _IEEE Transactions on Microwave Theory and Techniques_ , vol. 27, no. 12, pp. 987–993, 1979, doi: 10.1109/TMTT.1979.1129778.
  * [3] R. Marks, “A multiline method of network analyzer calibration,” _IEEE Transactions on Microwave Theory and Techniques_ , vol. 39, no. 7, pp. 1205–1215, 1991, doi: 10.1109/22.85388.
  * [4] H.-J. Eul and B. Schiek, “Thru-match-reflect: One result of a rigorous theory for de-embedding and network analyzer calibration,” in _1988 18th European Microwave Conference_ , 1988, doi: 10.1109/EUMA.1988.333924. pp. 909–914.
  * [5] L. Hayden, “An enhanced line-reflect-reflect-match calibration,” in _2006 67th ARFTG Conference_ , 2006, doi: 10.1109/ARFTG.2006.4734364. pp. 143–149.
  * [6] R. B. Marks, “Formulations of the basic vector network analyzer error model including switch-terms,” in _50th ARFTG Conference Digest_ , vol. 32, 1997, doi: 10.1109/ARFTG.1997.327265. pp. 115–126.
  * [7] S. Vandenberghe, D. Schreurs, G. Carchon, B. Nauwelaers, and W. De Raedt, “Identifying error-box parameters from the twelve-term vector network analyzer error model,” in _60th ARFTG Conference Digest, Fall 2002._ , 2002, doi: 10.1109/ARFTGF.2002.1218697. pp. 157–165.
  * [8] L. Hayden, “Vna error model conversion for n-port calibration comparison,” in _2007 69th ARFTG Conference_ , 2007, doi: 10.1109/ARFTG.2007.5456343. pp. 1–10.
  * [9] J. A. Jargon, D. F. Williams, and A. Sanders, “The relationship between switch-term-corrected scattering-parameters and wave-parameters measured with a two-port vector network analyzer,” _IEEE Microwave and Wireless Components Letters_ , vol. 28, no. 10, pp. 951–953, 2018, doi: 10.1109/LMWC.2018.2867076.
  * [10] J. A. Jargon and R. B. Marks, “Two-tier multiline trl for calibration of low-cost network analyzers,” in _46th ARFTG Conference Digest_ , vol. 28, 1995, doi: 10.1109/ARFTG.1995.327126. pp. 1–8.
  * [11] A. Rumiantsev and N. Ridler, “VNA calibration,” _IEEE Microwave Magazine_ , vol. 9, no. 3, pp. 86–99, jun 2008, doi: 10.1109/mmm.2008.919925.
  * [12] G. Strang, “The fundamental theorem of linear algebra,” _The American Mathematical Monthly_ , vol. 100, no. 9, pp. 848–855, nov 1993, doi: 10.2307/2324660.
  * [13] A. Arsenovic _et al._ , “scikit-rf: An open source python package for microwave network creation, analysis, and calibration [speaker’s corner],” _IEEE Microwave Magazine_ , vol. 23, no. 1, pp. 98–105, 2022, doi: 10.1109/MMM.2021.3117139.
  * [14] Z. Hatab, M. Gadringer, and W. Bösch, “Improving the reliability of the multiline trl calibration algorithm,” in _2022 98th ARFTG Microwave Measurement Conference (ARFTG)_ , 2022, doi: 10.1109/ARFTG52954.2022.9844064. pp. 1–5.