2022

[1]Hana Krásná

[1]Department of Geodesy and Geoinformation, Technische Universität Wien,
Wiedner Hauptstraße 8-10/E120.4, Vienna, 1040, Austria

2]United States Naval Observatory, USA

3]South African Radio Astronomy Observatory, South Africa

4]Jet Propulsion Laboratory, California Institute of Technology, USA

# The K-band (24 GHz) Celestial Reference Frame determined from Very Long
Baseline Interferometry sessions conducted over the past 20 years

<EMAIL_ADDRESS>David Gordon  Aletha de Witt  Christopher S. Jacobs
*  [  [  [

###### Abstract

The third realization of the International Celestial Reference Frame (ICRF3)
was adopted in August 2018 and includes positions of extragalactic objects at
three frequencies: 8.4 GHz, 24 GHz, and 32 GHz. In this paper, we present
celestial reference frames estimated from Very Long Baseline Interferometry
measurements at K-band (24 GHz) including data until June 2022. The data set
starts in May 2002 and currently consists of more than 120 24h observing
sessions performed over the past 20 years. Since the publication of ICRF3, the
additional observations of the sources during the last four years allow
maintenance of the celestial reference frame and more than 200 additional
radio sources ensure an expansion of the frame. A study of the presented
solutions is carried out helping us to understand systematic differences
between the astrometric catalogs and moving us towards a better next ICRF
solution. We compare K-band solutions (VIE-K-2022b and USNO-K-2022July05)
computed by two analysts with two independent software packages (VieVS and
Calc/Solve) and describe the differences in the solution strategy. We assess
the systematic differences using vector spherical harmonics and describe the
reasons for the most prominent ones.

###### keywords:

Very Long Baseline Interferometry, Celestial Reference Frame, K-band

##  1 Introduction

The current International Celestial Reference Frame (ICRF3; Charlot et al,
2020) is the third realization of the International Celestial Reference System
adopted by the International Astronomical Union in August 2018. The ICRF3 is
the first multi-wavelength radio frame since it contains positions of active
galactic nuclei (AGN) observed with Very Long Baseline Interferometry (VLBI)
at 2.3 and 8.4 GHz (S/X-band), 24 GHz (K-band), and 8.4 and 32 GHz (X/Ka-
band). The three components differ as shown by several statistical indicators
(e.g., data span, number of sources, coordinate uncertainty, error ellipse)
and each of them faces different challenges. In 2018 IAU Resolution B2, “On
The Third Realization of the International Celestial Reference Frame,” (ICRF3
working group, 2018) recommended that appropriate measures should be taken to
both maintain and improve ICRF3. In response, this paper concentrates on the
two main challenges in improving the accuracy of the celestial reference frame
observed at K-band (K-CRF) which are (1) observations at a single frequency
requiring an external ionospheric calibration and (2) the lack of a uniform
global terrestrial network causing a non-optimal observation geometry. Our
main goal is to assess systematic differences in the K-CRF solutions which are
computed at two VLBI analysis centers: at TU Wien with VLBI software package
VieVS (Böhm et al, 2018) and at the United States Naval Observatory (USNO)
with Calc/Solve. We also compare these two frames to the ICRF3 using vector
spherical harmonics (VSH) which provides information about systematic
differences between pairs of astrometric catalogs and we investigate the
possible reasons for the estimated differences.

##  2 Data and solution setup

###  2.1 Data description

Table 1: Overview of sessions included in our solutions listed with recording rate. time span | session code | data rate [Mbps]
---|---|---
| Northern (VLBA) sessions |
05/2002 - 12/2008 | BR079a-c, BL115a-c, BL122a-d, BL151a-b | 128
06/2006 - 10/2006 | BP125a-c | 256
12/2015 - 10/2019 | BJ083a-d, UD001a-x, UD009a-o | 2048
11/2019 - 06/2022 | UD009p-z, UD009aa-ah, UD015a-l | 4096
| Southern sessions |
05/2014 - 07/2016 | KS1401, KS1601 | 1024
11/2016 - 02/2021 | KS1603, KS1702-KS2102 | 2048
|  |

Table 2: Selected models and parametrization in VIE-K-2022b and
USNO-K-2022July05. Values in parentheses represent the applied constraints.
The abbreviation pwlo stands for piecewise linear offset.

A priori modeling | VIE-K-2022b | USNO-K-2022July05
---|---|---
ionosphere maps | CODE time series (Schaer, 1999) | 2 h average JPL maps
ionospheric mf coefficients | MSLM, $k$ = 1, $\Delta H$ = 56.7 km, $\alpha$ = 0.9782 | 2-D thin shell, MSLM
hydrostatic delay | in situ pressure (Saastamoinen, 1972) | in situ pressure (Saastamoinen, 1972)
hydrostatic + wet mf | VMF3 (Landskron and Böhm, 2018) | VMF1 (Böhm et al, 2006)
hydrostatic gradients | DAO (MacMillan and Ma, 1997) | DAO (MacMillan and Ma, 1997)
precession/nutation model | IAU 2006/2000A | IAU 2006/2000A
celestial pole offsets (CPO) | IERS Bulletin A, https://maia.usno.navy.mil/ser7/finals2000A.all | none
Parametrization |  |
zenith wet delay | 30 min pwlo (1.5 cm/30 min) | 30 min pwlo (1.5 cm/h)
tropo. grad.: VLBA | 3 h pwlo (0.5 mm/3 h) | 6 h pwlo (0.5mm, 2 mm/day)
tropo. grad.: KS | fixed to a priori | fixed to a priori
CPO: VLBA | 24 h pwlo (0.1 µas/24 h) | offset at midpoint of the session
CPO: KS | fixed to a priori | fixed to a priori
weighting | elevation-dependent (Gipson et al, 2008) | baseline-dependent
|  |

The celestial reference frames introduced in this paper are computed from
$1.96\cdot 10^{6}$ group delays observed at K-band in the VLBI sessions listed
in Table 1. This data set was acquired mainly with the Very Long Baseline
Array (VLBA) starting in May 2002. The first sessions belong to programs
carried out by Lanyi et al (2010) and Petrov et al (2011). All sessions up to
May 2018 are part of the current ICRF at K-band, ICRF3-K. The VLBA (Napier,
1995), because its sites are limited to U.S. territory, does not allow
observations of sources with declinations below $-46^{\circ}$. Therefore,
southern K-band sessions (KS) were organized starting in May 2014. The vast
majority of southern observations are from single baseline sessions between
the HartRAO 26m (South Africa) and the Hobart 26m (Tasmania, Australia) with
the exception of one session involving the Tianma 65m (near Shanghai, China)
and four sessions augmented with the Tidbinbilla 70m telescope (near Canberra,
Australia). Of all the sources, 913 were observed in VLBA sessions, 328 were
observed in southern hemisphere sessions, and 206 were observed in both types
between $-46^{\circ}$ and $+39^{\circ}$ declination.

In Fig. 1 we show the number of observations conducted after the ICRF3 K-band
data cutoff on 5 May 2018 until June 2022 divided into three groups: (a)
observations to ICRF3-K defining sources, (b) observations to ICRF3-K non-
defining sources and (c) observations to sources which are not included in
ICRF3-K. The consequence of using mainly the VLBA network for the K-band
observations is the lack of observations of the deep south sources which is
currently amplified by the technical problems of Hobart26 since March 2021.
The low number of new observations (under 100) of the deep south sources since
the ICRF3 release is seen in all three plots of Fig. 1.

Figure 1: Number of observations after the ICRF3-K data cutoff in May 2018
until June 2022. The sources are divided into three groups: ICRF3-K defining
sources (top), ICRF3-K non-defining sources (middle), and sources not included
in ICRF3-K (bottom).

###  2.2 Setup of solutions

The treatment of the K-band VLBI observations in the VieVS solution
(VIE-K-2022b) is similar to the S/X solution VIE2022b computed at the VIE
Analysis Center111https://www.vlbi.at of the International VLBI Service for
Geodesy & Astrometry. A detailed description of the setup and applied
theoretical models during the analysis are given in Krásná et al (2022). In
Table 2 we highlight models used in VIE-K-2022b and the USNO Calc/Solve
solution USNO-K-2022July05222latest version at
https://crf.usno.navy.mil/data_products/RORFD/Quarterly/current//USNO_Kband_source_positions.iers
relevant to the presented investigations. While S/X frames calibrate the
ionosphere directly from their dual-band data, K-band ionospheric effects
require external calibration data. Specifically, K-band systems at the VLBA
and the southern stations currently lack the complementary lower band needed
for a dual-band ionospheric calibration, therefore the frequency-dependent
delay coming from the dispersive part of the atmosphere has to be described by
external models. In both K-band solutions presented here, ionospheric maps
derived from Global Navigation Satellite System (GNSS) are applied. In
VIE-K-2022b, global ionospheric maps provided by the Center for Orbit
Determination in Europe (CODE; Schaer, 1999)333http://ftp.aiub.unibe.ch/CODE/
are used with a time spacing of two hours from 05/2002 until 05/2014, and of
one hour since that date. In USNO-K-2022July05, global ionospheric maps
computed at the Jet Propulsion Laboratory (JPL) with two hours resolution are
applied.
The alignment of the Terrestrial Reference Frame (TRF) is done by applying the
No-Net-Translation (NNT) and No-Net-Rotation (NNR) conditions to the station
position and velocity parameters in the global normal matrix. In VIE-K-2022b,
the conditions are applied to all VLBA telescopes but one (MK-VLBA) with
respect to the ITRF2020. In USNO-K-2022July05, the NNT/NNR condition is used
w.r.t. a TRF solution based on ITRF2014 applied to all participating antennas
except MK-VLBA (position discontinuity due to an Earthquake on June 15, 2006)
and TIDBIN64 (limited number of observations).

The common practice for the rotational alignment of a new celestial reference
frame to the current official one is to apply a three-dimensional constraint
to the defining sources. In both solutions, ICRF3-SX is used as a priori
celestial reference frame and the galactic acceleration correction is modeled
with the adopted ICRF3 value of 5.8 µas/yr for the amplitude of the solar
system barycenter acceleration vector for the epoch 2015.0. Datum definition
of the CRFs is accomplished by the unweighted NNR (Jacobs et al, 2010) w.r.t.
287 (VIE-K-2022b) and 258 (USNO-K-2022July05) defining ICRF3-SX sources.

##  3 Results

We analyze the estimated VIE-K-2022b and USNO-K-2022July05 frames in terms of
the vector spherical harmonics decomposition (VSH; Mignard and Klioner, 2012;
Titov, O. and Lambert, S., 2013; Mayer and Böhm, 2020) w.r.t. ICRF3-SX which
allows studying possible systematic differences between the catalogs. Prior to
the comparison, outliers—defined as AGN with an angular separation greater
than 5 mas from their ICRF3-SX position—were removed. In both solutions, there
are four outlier sources: 0134+329 (3C48), 0316+162 (CTA21), 0429+415 (3C119),
and 2018+295. Note that large position changes for 3C48 and CTA21 were found
at X-band in observations made after the ICRF3 release and are reported by
Frey and Titov (2021) and Titov et al (2022). The number of remaining common
sources is 993 in VIE-K-2022b and 995 in USNO-K-2022July05. The two sources
(0227-542 and 0517-726) missing in VIE-K-2022b have 3 and 4 observations in
USNO-K-2022July05. In VIE-K-2022b these observations were removed based on an
outlier check of individual observations during the single session analysis.

The VSH are obtained with a least squares adjustment where the weight matrix
contains inflated formal errors of the source coordinates. Similar to ICRF3-K,
the formal errors of the source coordinates in both catalogs are inflated by a
factor of 1.5, and a noise floor of 30 and 50 µas in quadrature is added to
right ascension and declination, respectively. Table 3 summarizes the first
order and second degree and order VSH, i.e., rotation ($R_{1},R_{2},R_{3}$),
dipole ($D_{1},D_{2},D_{3}$), and ten coefficients ($a$) for the quadrupole
harmonics of magnetic ($m$) and electric ($e$) type. All three rotation angles
between the VIE-K-2022b and ICRF3-SX axes are within their formal errors and
the angles do not exceed 8 µas. The largest angle ($16\pm 10$ µas) between
USNO-K-2022July05 and ICRF-SX is around the y-axis ($R_{2}$). The selection of
defining sources for the NNR constraint influences the mutual rotations of two
catalogs (cf. Section 3.1 for more details). The three dipole parameters
represent the distortion as a flow from a source to a sink located at two
opposite poles. The $D_{3}$ term ($-4\pm 10$ µas in VIE-K-2022b and $60\pm 9$
µas in USNO-K-2022July05) is susceptible to imperfect modeling of equatorial
bulges in the ionospheric and tropospheric calibrations (cf. Section 3.2). The
zonal quadrupole terms $a_{2,0}^{e}$ and $a_{2,0}^{m}$ reflect north-south
asymmetries. Their values w.r.t. ICRF3-SX reach $-3\pm 12$ µas and $-36\pm 7$
µas in VIE-K-2022b, and $-46\pm 11$ µas and $1\pm 7$ µas in USNO-K-2022July05,
respectively (cf. Section 3.3).

Table 3: VSH parameters up to degree and order two for VIE-K-2022b and USNO-K-2022July05 w.r.t. ICRF3-SX (after eliminating four outliers from the solutions). [µas] | VIE-K-2022b | USNO-K-2022July05
---|---|---
$R_{1}$ | $-1~{}\pm 10$ | $-4~{}\pm 10$
$R_{2}$ | $-8~{}\pm 10$ | $-16~{}\pm 10$
$R_{3}$ | $+0~{}\pm~{}6$ | $-11~{}\pm~{}6$
$D_{1}$ | $-17~{}\pm~{}9$ | $-5~{}\pm~{}9$
$D_{2}$ | $-15~{}\pm~{}9$ | $+9~{}\pm~{}9$
$D_{3}$ | $-4~{}\pm 10$ | $+60~{}\pm~{}9$
$a_{2,0}^{e}$ | $-3~{}\pm 12$ | $-46~{}\pm 11$
$a_{2,0}^{m}$ | $-36~{}\pm~{}7$ | $+1~{}\pm~{}7$
$a_{2,1}^{e,Re}$ | $-19~{}\pm 10$ | $-13~{}\pm 10$
$a_{2,1}^{e,Im}$ | $-21~{}\pm 11$ | $-26~{}\pm 11$
$a_{2,1}^{m,Re}$ | $-13~{}\pm 11$ | $+13~{}\pm 10$
$a_{2,1}^{m,Im}$ | $-12~{}\pm 11$ | $-6~{}\pm 11$
$a_{2,2}^{e,Re}$ | $+1~{}\pm~{}4$ | $+3~{}\pm~{}4$
$a_{2,2}^{e,Im}$ | $+8~{}\pm~{}4$ | $+3~{}\pm~{}4$
$a_{2,2}^{m,Re}$ | $+12~{}\pm~{}5$ | $+23~{}\pm~{}5$
$a_{2,2}^{m,Im}$ | $+6~{}\pm~{}5$ | $+4~{}\pm~{}5$

###  3.1 Defining sources

During the development of the ICRF3 a new set of sources observed at S/X-band
was selected for defining the rotational alignment. This set of defining
sources was based on three selection criteria in order to align the S/X-frame
with its predecessor, the ICRF2 (Fey et al, 2015). These criteria were: (1)
the overall sky distribution of the defining sources, (2) the position
stability of the individual sources, and (3) the compactness of their
structures (Charlot et al, 2020). For the alignment of the K-band reference
frame ICRF3-K, a subset of 193 sources out of the set of 303 ICRF3-S/X
defining sources—based mainly on the number of available K-band
observations—was used. In Fig. 2 we show the distribution of the ICRF3-SX
defining sources and highlight the ICRF3-K defining sources with yellow color.
In the solutions VIE-K-2022b (red crosses) and USNO-K-2022July05 (blue dots)
we take advantage of the additional observations gained after the ICRF3
release and choose the defining sources independently of the ICRF3-K ones. The
current analysis of available sessions shows that there are no K-band
observations of four ICRF3-S/X defining sources: 0044-846, 0855-716, 1448-648,
1935-692. This means, that 299 out of the 303 ICRF3-SX defining sources are
observed in K-band (considering June 2022 to be the cutoff date for K-band
observations). In VIE-K-2022b and USNO-K-2022July05 we apply different
strategies for the selection of defining sources.

Figure 2: Defining sources. The circles denote the 303 ICRF3-SX defining
sources. The subgroup of 193 yellow circles depicts the ICRF3-K defining
sources. Defining sources in VIE-K-2022b and USNO-K-2022July05 are red crosses
and blue dots, respectively.

At TU Wien, we first computed a K-CRF solution from VLBA sessions only. We
found 12 AGN (0038-326, 0227-369, 0316-444, 0437-454, 0743-006, 1143-245,
1606-398, 1929-457, 1937-101, 2036-034, 2111+400, 2325-150) among the 303
ICRF3-SX defining sources whose angular separation in this VLBA-only K-CRF
solution is greater than 0.5 mas from their ICRF3-SX position and those are
dropped from the NNR condition in VIE-K-2022b. All ICRF3-SX defining sources
observed in the KS sessions only are kept in the NNR in VIE-K-2022b.

In USNO-K-2022July05 the following sources were excluded from the defining
set: 0700-465, 0742-562, 0809-493 and 0918-534 since they show offsets of
$0.5-1.5$ mas from their ICRF3-SX positions in recent USNO S/X solutions. An
additional 41 sources, mostly in the deep south, were also excluded from the
NNR condition because they had either very few or no observations.

The rotation angles in Table 3 show that the incorporation of the deep south
sources in the alignment condition makes the adjustment more robust and keeps
the estimated K-CRF solution slightly closer to the a priori one.

###  3.2 Ionospheric mapping function

The global ionosphere maps provide the Vertical Total Electron Content (VTEC).
The conversion from VTEC to the Slant Total Electron Content (STEC) at an
elevation angle ($\epsilon$) of the VLBI observations at the telescope is done
by the ionospheric mapping function (mf, $M$). In VIE-K-2022b we apply the
thin shell ionospheric mf introduced by Schaer (1999) and recently discussed
in detail by Petrov (2023):

$M(\epsilon)=k\cdot\frac{1}{\sqrt{1-\Big{(}\frac{R_{E}}{R_{E}+H_{i}+\Delta
H}\Big{)}^{2}\cdot\cos^{2}\alpha\epsilon}},$ (1)

where $k$ is a scaling factor, $R_{E}$ = 6371 km stands for the Earth’s base
radius, $H_{i}$ = 450 km is the height of the spherical single layer, $\Delta
H$ represents an increment in the ionosphere height, and $\alpha$ is a
correction factor to the elevation angle. In the default VIE-K-2022b solution
we apply: $k$ = 1, $\Delta H$ = 56.7 km, $\alpha$ = 0.9782 which is denoted as
Modified Single-Layer Model
(MSLM)444http://ftp.aiub.unibe.ch/users/schaer/igsiono/doc/mslm.pdf mapping
function and claimed to be the best fit with respect to the JPL extended slab
model mapping function. This parameter setting is recommended e.g. by Feltens
et al (2018), Wielgosz et al (2018), and references therein. The standard
Single Layer Model (SLM) mapping function is achieved with the parameters: $k$
= 1, $\Delta H$ = 0 km, and $\alpha$ = 1. Following the discussion in Petrov
(2023), we calculated two more solutions with different ionospheric mf
parametrizations based on MSLM with different values of $\Delta H$ and $k$
(i.e., iono3 and iono4) as summarized in Table 4.

In order to quantify the effect of the modified ionospheric mapping function
on the K-CRF solution, we calculated VSH for each solution w.r.t. ICRF-SX
(Fig. 3). Changes in the three mf parameters ($k$, $\Delta H$, $\alpha$)
influence the terms $D_{3}$ and $a_{2,0}^{e}$, which are sensitive to the
equatorial bulge and north-south asymmetries, as mentioned earlier. The best
fit to the ICRF3-SX is achieved with the MSLM mapping function applied in
VIE-K-2022b where these parameters are negligibly small ($-4\pm 10$ µas and
$-3\pm 12$ µas, respectively). On the other hand, in iono4 (where a scale
factor $k$ = 0.85 is applied to MSLM), the difference w.r.t. ICRF3-SX in
$D_{3}$ and $a_{2,0}^{e}$ increases to $42\pm 10$ µas and $-15\pm 12$ µas,
respectively. In Fig. 4 we plot the differences in declination between the
four discussed solutions w.r.t. ICRF3-SX over declination for individual
sources. The smoothed curves are computed as moving averages with a Gaussian
kernel and plotted with color coding identical to Fig. 3. The positive
systematic difference in the declination estimates w.r.t. ICRF3-SX, appearing
approximately between -40∘ and -10∘ declination, reaches its maximum of 63 µas
for -26∘ declination in VIE-K-2022b with applied MSLM mapping function (blue
curve).

Table 4: Parameters of the ionospheric mapping function and the resulting VSH parameters $D_{3}$ and $a_{2,0}^{e}$.  | $k$ | $\Delta H$ | $\alpha$ | $D_{3}$ | $a_{2,0}^{e}$
---|---|---|---|---|---
| [-] | [km] | [-] | [µas] | [µas]
MSLM | 1 | 56.7 | 0.9782 | $-4\pm 10$ | $-3\pm 12$
SLM | 1 | 0 | 1 | $-17\pm 10$ | $+2\pm 12$
iono3 | 1 | 150.0 | 0.9782 | $15\pm 10$ | $-10\pm 12$
iono4 | 0.85 | 56.7 | 0.9782 | $42\pm 10$ | $-15\pm 12$
Figure 3: Vector spherical harmonics of K-CRF solutions computed with
ionospheric mf MSLM (blue, VIE-K-2022b), SLM (red), iono3 (yellow), and iono4
(green) w.r.t. ICRF3-SX. Figure 4: Smoothed differences in declination from
K-CRF solutions computed with ionospheric mf MSLM (blue, VIE-K-2022b), SLM
(red), iono3 (yellow), and iono4 (green) w.r.t. ICRF3-SX. The black dots are
differences in the declination of individual sources in VIE-K-2022b w.r.t.
ICRF3-SX and their formal errors (in grey).

###  3.3 Systematic in elevation angles

Along with the ionospheric effects, the K-CRF suffers from an asymmetric
observing network geometry with 99% of the data being from the all-northern
VLBA. In Fig. 5 the percentage of observations from southern KS sessions for
individual sources in VIE-K-2022b is shown. The logarithmic color scale
highlights the fact that the number of observations to the sources with
declination higher than $-45^{\circ}$ builds only a tiny fraction of the total
number of observations, although the mutual sky visibility between southern
antennas HartRAO and Hobart allows observing sources up to approximately
$30^{\circ}$ declination. The mean percentage of observations from KS sessions
for sources with declination between $-15^{\circ}$ and $-45^{\circ}$ (area
with mainly yellow and green colors in Fig. 5) is 0.96%. The total number of
K-CRF observations (blue color) and the number of observations from the KS
(brown color) during individual years is plotted in Fig. 6. The numbers above
the columns give the percentage of observations from KS w.r.t. the total
number of observations within the individual year.

Figure 5: Percentage of observations from KS sessions among the total number
of K-CRF observations for individual sources in VIE-K-2022b. Figure 6: The
yearly distribution of the K-CRF observations used in our solutions. The
numbers above the columns show the percentage of observations from southern
sessions (KS, brown columns) w.r.t. the total number of observations (blue
columns) during the individual year.

In order to explore the resultant elevation-dependent effects, we characterize
the distribution of elevation angles at which the sources were observed. These
distributions vary due to both the geometry of the VLBA network and the fact
that we observe each source over a range of hour angles. First, we define a
parameter called airmass in order to quantify the approximate total pathlength
through the troposphere for each source—with the maximum at lower elevation
angles. It is computed for each observation from the whole data set with the
simplifying assumption of a flat slab atmosphere (ignoring the curvature of
the atmosphere over a spherical Earth):

$\mathrm{airmass}=\frac{1}{\sin(\epsilon_{1})}+\frac{1}{\sin(\epsilon_{2})},$
(2)

where $\epsilon$ is the elevation angle of the source at telescopes 1 and 2 of
the baseline. Next, we compute the median value over the individual
observations for each source and plot it with respect to the declination (Fig.
7) with the errors (in grey) obtained as standard deviations computed over the
individual airmass values for the particular source. The systematic increase
of the airmass parameter from 0∘ to $-45^{\circ}$ declination can lead to an
overestimation of the optimal data weights for VLBA observations in this
declination range when the larger noise of observations conducted at low
elevation angles is not considered.

Figure 7: Median airmass for individual sources computed over all their
observations in VIE-K-2022b.

To partly account for the overweighting of the low elevation scans (which
observe low declination sources in the mentioned area), elevation-dependent
weighting (Eq. (3); Gipson et al, 2008) in VIE-K-2022b is applied. In the
diagonal covariance matrix the measurement noise $\sigma^{2}_{m}$ is increased
by the squared elevation-dependent noise terms for telescopes 1 and 2:

$\sigma^{2}_{obs}=\sigma^{2}_{m}+\bigg{(}\frac{6~{}\texttt{ps}}{\sin(\varepsilon_{1})}\bigg{)}^{2}+\bigg{(}\frac{6~{}\texttt{ps}}{\sin(\varepsilon_{2})}\bigg{)}^{2}.$
(3)

Hence, sources between 0∘ to $-45^{\circ}$ declination obtain a lower weight
in the least squares adjustment and the resulting distortion of the celestial
reference frame is damped. For example, an observation conducted with two VLBA
antennas at the elevation angles of 15∘ has an airmass value of 8 (Eq. (2))
which corresponds in our data set to a source with a declination of about
$-40^{\circ}$ (Fig. 7). The additional noise added to the $\sigma^{2}_{m}$ of
this observation in quadrature is 33 ps (Eq. (3)) which decreases its weight
in the solution.

##  4 Conclusion

Recent K-CRF solutions computed at TU Wien
(VIE-K-2022b)555https://vlbi.at/data/analysis/ggrf/crf_vie2022b_k.txt and USNO
(USNO-K-2022July05) from single-frequency band VLBI observations (24 GHz)
until June 2022 were assessed. The vector spherical harmonics were computed
w.r.t. ICRF3-SX after eliminating four AGN as outliers. In VIE-K-2022b, all
rotation values are lower than 8 µas and have significance at the level of
their formal errors or less. With a single exception, all dipole and
quadrupole terms are within 20 µas with a marginal significance of two times
the formal error as maximum. The only quadrupole term above this limit is
$a_{2,0}^{m}$ ($-36\pm 7$ µas). We discussed two major challenges which limit
the accuracy of the current K-band VLBI solutions: external ionospheric
corrections and the non-uniform observing network geometry—especially the lack
of observations in the deep south. We show that the choice of ionospheric
mapping function parameters influences the dipole, $D_{3}$, and quadrupole
terms $a_{2,0}^{e}$. Because 99% of the data is observed with the all-northern
VLBA sources between 0∘ and $-45^{\circ}$ declination have a monotonic
decrease in median elevation angle of observation making our solution
vulnerable to atmospheric mis-modeling. We reduced sensitivity of the
VIE-K-2022b solution to the effect of this observing geometry bias by
computing elevation-dependent weighting to downweight low elevation
observations. Future work will focus on improving the geometry of the K-band
observing network, improving the modeling of atmospheric effects, and
improving solution weighting schemes.

##  5 Declarations

Ethics approval and consent to participate Not applicable. Consent for
publication Not applicable. Competing interests There are no relevant
financial or non-financial competing interests to report. Funding We
acknowledge our respective sponsors: SARAO/HartRAO is a facility of the
National Research Foundation (NRF) of South Africa. Portions of this work were
done at the Jet Propulsion Laboratory, California Institute of Technology
under contract with NASA (contract no. 80NM0018D0004). Portions of this work
were sponsored by the Radio Optical Reference Frame Division of the U.S. Naval
Observatory. This work supports USNO’s ongoing research into the celestial
reference frame and geodesy. Authors’ contributions HK wrote the manuscript,
analyzed the VLBI data and created the VIE-K solutions. DG prepared the vgosDB
databases and computed the USNO-K solution. AdW is the PI of the VLBI K-band
group and the leader by planning of the VLBI K-band observations. CJ proposed
the concept of the paper and contributed to the analysis of data. All authors
contributed to regular discussions and interpretations of results. They read
and commented the final paper. Acknowledgments The authors appreciate comments
provided by three anonymous reviewers. HK thanks Leonid Petrov (NASA GSFC) for
fruitful discussions about single band astrometry. The authors gratefully
acknowledge the use of the VLBA under the USNO’s time allocation.

## References

  *   * Böhm et al (2006) Böhm J, Werl B, Schuh H (2006) Troposphere mapping functions for GPS and very long baseline interferometry from European Centre for Medium-Range Weather Forecasts operational analysis data. J Geophys Res: Solid Earth 111(B2). 10.1029/2005JB003629
  * Böhm et al (2018) Böhm J, Böhm S, Boisits J, et al (2018) Vienna VLBI and Satellite Software (VieVS) for Geodesy and Astrometry. Publ Astron Soc Pac 130(986):044,503. 10.1088/1538-3873/aaa22b
  * Charlot et al (2020) Charlot P, Jacobs CS, Gordon D, et al (2020) The third realization of the International Celestial Reference Frame by very long baseline interferometry. Astron Astrophys 644:A159. 10.1051/0004-6361/202038368
  * Feltens et al (2018) Feltens J, Bellei G, Springer T, et al (2018) Tropospheric and ionospheric media calibrations based on global navigation satellite system observation data. J Space Weather Space Clim 8:A30. 10.1051/swsc/2018016
  * Fey et al (2015) Fey AL, Gordon D, Jacobs CS, et al (2015) The second realization of the International Celestial Reference Frame by Very Long Baseline Interferometry. Astron J 150(2):58. 10.1088/0004-6256/150/2/58
  * Frey and Titov (2021) Frey S, Titov O (2021) Change in the radio structure and position of the quasar CTA 21. Research Notes of the AAS 5(3):60. 10.3847/2515-5172/abf123
  * Gipson et al (2008) Gipson J, MacMillan D, Petrov L (2008) Improved Estimation in VLBI through Better Modeling and Analysis. In: Finkelstein A, Behrend D (eds) IVS GM Proceedings, pp 157–162, URL https://ui.adsabs.harvard.edu/abs/2008mefu.conf..157G/abstract
  * ICRF3 working group (2018) ICRF3 working group (2018) IAU Resolution B2 on The Third Realization of the International Celestial Reference Frame. In: Proceedings of the XXX IAU General Assembly. IAU, Vienna, Austria, URL https://www.iau.org/static/resolutions/IAU2018_ResolB2_English.pdf
  * Jacobs et al (2010) Jacobs C, Heflin M, Lanyi G, et al (2010) Rotational alignment altered by source position correlations. In: Behrend D, Baver KD (eds) IVS GM Proceedings. NASA, pp 305–309, URL https://ivscc.gsfc.nasa.gov/publications/gm2010/jacobs2.pdf
  * Krásná et al (2022) Krásná H, Baldreich L, Böhm J, et al (2022) VLBI celestial and terrestrial reference frames VIE2022b 10.48550/arXiv.2211.07338
  * Landskron and Böhm (2018) Landskron D, Böhm J (2018) VMF3/GPT3: refined discrete and empirical troposphere mapping functions. J Geodesy 92(4):349––360. 10.1007/s00190-017-1066-2
  * Lanyi et al (2010) Lanyi GE, Boboltz DA, Charlot P, et al (2010) The celestial reference frame at 24 and 43 GHz. I. Astrometry. Astron J 139(5):1695. 10.1088/0004-6256/139/5/1695
  * MacMillan and Ma (1997) MacMillan DS, Ma C (1997) Atmospheric gradients and the VLBI terrestrial and celestial reference frames. Geophys Res Lett 24(4):453–456. 10.1029/97GL00143
  * Mayer and Böhm (2020) Mayer D, Böhm J (2020) Comparing vienna crf solutions to gaia-crf2. In: Freymueller JT, Sánchez L (eds) Beyond 100: The Next Century in Geodesy. Springer International Publishing, Cham, pp 21–28, 10.1007/1345_2020_99
  * Mignard and Klioner (2012) Mignard F, Klioner S (2012) Analysis of astrometric catalogues with vector spherical harmonics. Astron Astrophys 547:A59. 10.1051/0004-6361/201219927
  * Napier (1995) Napier P (1995) VLBA Design. In: Zensus J, Diamond P, Napier P (eds) Very Long Baseline Interferometry and the VLBA, in ASP Conference Series, vol 82. Astronomical Society of the Pacific, p 59–72, URL https://articles.adsabs.harvard.edu/pdf/1995ASPC...82...59N
  * Petrov (2023) Petrov L (2023) Single-band VLBI Absolute Astrometry. Astron J 165(4):183. 10.3847/1538-3881/acc174
  * Petrov et al (2011) Petrov L, Kovalev YY, Fomalont EB, et al (2011) The Very Long Baseline Array Galactic Plane Survey—VGaPS. Astron J 142(2):35. 10.1088/0004-6256/142/2/35
  * Saastamoinen (1972) Saastamoinen J (1972) Introduction to practical computation of astronomical refraction. Bull Geodesique 106:383–397. 10.1007/BF02522047
  * Schaer (1999) Schaer S (1999) Mapping and predicting the Earth’s ionosphere using the global positioning system. PhD thesis, University of Berne
  * Titov et al (2022) Titov O, Frey S, Melnikov A, et al (2022) Unprecedented change in the position of four radio sources. Mon Not R Astron Soc 512(1):874–883. 10.1093/mnras/stac038
  * Titov, O. and Lambert, S. (2013) Titov, O., Lambert, S. (2013) Improved VLBI measurement of the solar system acceleration. Astron Astrophys 559:A95. 10.1051/0004-6361/201321806
  * Wielgosz et al (2018) Wielgosz P, Milanowska B, Krypiak-Gregorczyk A, et al (2018) Validation of gnss-derived global ionosphere maps for different solar activity levels: case studies for years 2014 and 2018. GPS Solutions 25(103). 10.1007/s10291-021-01142-x


#  Temporal and Spatial Online Integrated Calibration for Camera and LiDAR
††thanks: This work was supported by the National High Technology Research and
Development Program of China under Grant No. 2018YFE0204300, and the National
Natural Science Foundation of China under Grant No. U1964203, and sponsored by
Tsinghua University-Didi Joint Research Center for Future Mobility.

Shouan Wang1, Xinyu Zhang1,∗, GuiPeng Zhang1, Yijin Xiong1, Ganglin Tian1,
Shichun Guo1, Jun Li1, Pingping Lu2, Junqing Wei3 and Lei Tian4 1The authors
are with the State Key Laboratory of Automotive Safety and Energy, the School
of Vehicle and Mobility, Tsinghua University, 100084.2Pingping Lu is with
University of Michigan.3Junqing Wei is with DiDi Co. Ltd.4Lei Tian is with
SINOTRUK Co. Ltd.*Author to whom correspondence should be addressed. (e-mail:
<EMAIL_ADDRESS>
###### Abstract

While camera and LiDAR are widely used in most of the assisted and autonomous
driving systems, only a few works have been proposed to associate the temporal
synchronization and extrinsic calibration for camera and LiDAR which are
dedicated to online sensors data fusion. The temporal and spatial calibration
technologies are facing the challenges of lack of relevance and real-time. In
this paper, we introduce the pose estimation model and environmental robust
line features extraction to improve the relevance of data fusion and instant
online ability of correction. Dynamic targets eliminating aims to seek optimal
policy considering the correspondence of point cloud matching between adjacent
moments. The searching optimization process aims to provide accurate
parameters with both computation accuracy and efficiency. To demonstrate the
benefits of this method, we evaluate it on the KITTI benchmark with ground
truth value. In online experiments, our approach improves the accuracy by
38.5% than the soft synchronization method in temporal calibration. While in
spatial calibration, our approach automatically corrects disturbance errors
within 0.4 second and achieves an accuracy of 0.3-degree. This work can
promote the research and application of sensor fusion.

###### Index Terms:

Calibration, pose estimation, dynamic targets elimination, autonomous driving

##  I INTRODUCTION

Figure 1: Line features are extracted from a point cloud (a) at moment
$t+e_{t}$ and an image (b) at moment $t$ of the same scene. Through pose
estimation model and calculation of correspondence between 2D and 3D line
features, the projected point cloud line features (c) are well aligned with
the image line features.

There have been many studies on multi-sensor fusion in autonomous vehicles in
recent years, benefiting from the multi-sensors equipped in autonomous
vehicles. Sensor fusion is an essential and critical strategy to enhance the
performance and ensure the reliability of perception modules in ADAS, such as
in target classification[1], segmentation[2], SLAM[3], etc. Both temporal and
spatial calibrations are the most critical procedures for sensor fusion. There
have been a number of methods and systems for temporal calibrations and
spatial calibrations.

For temporal calibration, hardware-based approaches are widely adopted for its
accuracy. However, a unified clock signal and global timestamp are required
and sensor hardware trigger interface must be reserved, which is not
universally applicable to all devices; software-based approaches such as soft-
synchronous interpolation[4], values are aligned by curve fitting, which
brings accuracy problems because the data are theoretical values obtained by
fitting to the theoretical value.

Figure 2: The architecture of the proposed calibration system.

For spatial calibration, specially designed objects are required for
traditional manual calibration approaches, such as manually selected
points[5], which brings redundant manual operations. Long-time operation and
varying loads can lead to slight drifts and deviations in extrinsic
parameters. Current artificially designed targets[6] are utilized to calibrate
the extrinsic parameters in automatic calibration works, but it’s only at the
stage of the laboratory. Some feature-based calibration method[7] utilizes
edge features to compute the extrinsic parameters. However, these features are
not well corresponded to each other under some scenarios.

The current temporal and spatial calibration modules are not well integrated
and associated, and only a few works were proposed to solve the temporal and
spatial calibration integration for camera and LiDAR through calibration
boards while ignoring the online performance.

Therefore, in this paper we aim to combine the temporal-spatial calibration
into online integrated calibration. For this purpose, sensors’ pose estimation
model between adjacent moments is introduced for online integrated
calibration, and line features are used to correct the error of the pose
estimation model.

The main contributions of this paper are as follows:

1) A novel pose estimation model is proposed to decrease time delay and
extrinsic offset in temporal and spatial online integrated calibration
procedures.

2) We introduce a method for eliminating dynamic point clouds which only used
the association of adjacent two frames instead of detection boxes based on
prior information

3) We introduced a searching optimization method to improve the optimization
efficiency and proposed a new set of accuracy assessment metrics to evaluate
the accuracy of the calibration results.

##  II RELATED WORKS

In this section, we discuss the development of temporal and spatial
calibrations. Temporal calibration and spatial calibration are poorly coupled
and generally separated. For the temporal calibration part, hardware-based
methods have been widely applied, due to their exceptional accuracy. The
proposal of TriggerSync[8] represented the high accuracy of the hardware
triggered method, however, it was low robust to the false correlation of
trigger pulses due to unexpected delays, and more interfaces were required. To
solve the problem of an unexpected delay, the Global Position System(GPS)[9]
signal was proposed as the timing signal. However, in the GPS rejected region,
timing accuracy was seriously affected by the signal problem. To overcome this
drawback, Marsel Faizullin et al.[10] proposed a microcontroller-based
platform for GPS signal emulation that can be extended to arbitrary multi-
sensors beyond LiDAR-imu. Hannes Sommer et al.[11] used simple external
devices, including LED (for camera) or photodiodes (for LiDAR), to achieve
temporal calibration. This was done at the cost of additional equipment in
exchange for a high rate and high precision estimate of time offset. Hardware-
triggered interfaces or support from other hardware were required by all of
the above approaches, which lead to their limited use in all scenarios.
Instead, software-based methods were gradually popular because of its low
hardware requirements. The message_filter method was based on Robot Operating
System(ROS), however, it only adopted a simple strategy to match timestamps
when dealing with sensor data of different frequencies, which performed poorly
in high-precision data fusion scenarios. bidirectional synchronization methods
such as TICSync[12], PTP or NTP, performed well but require sensor firmware
support, which was rare and expensive.

For the spatial calibration part, checkerboard and selected points were first
chosen to solve the calibration problem between LiDAR and camera[13]. The main
limitation of the above approaches was manual and time-consuming. To overcome
the above limitations, several ways were proposed by some researchers to
calibrate the extrinsic parameters more intelligently. Special artificial
calibration targets, such as spheres [14], were introduced to avoid manually
selecting. However, these methods still shared the limitation of specially
designed targets and man-made target errors. To solve this problem, targetless
methods were proposed. [15] was based on the maximization of mutual
information obtained between surface intensities measured. In [16], they
presented automatic targetless calibration methods based on hand-eye
calibration, which only required three-dimensional point clouds and camera
images to compute motion information, with a rough calibration result.

In contrast to the above separate approaches for temporal and spatial
calibration, some studies focused on integrated calibration. Chia-Le Lee et
al.[17] proposed a method to achieve temporal and extrinsic calibration
between radar and LiDAR based on a specified triangular plate target. And Z.
Taylor et al.[18] proposed a targetless self-calibration method that relied on
the same motion observed in different sensors to find their relative offsets.
All these methods prompted us to wonder whether there is an approach to solve
the above problems by online integrated temporal and spatial calibration based
on environmental features. Therefore, we propose a method based on the line-
based extrinsic calibration method [19], extending the applicability of the
extrinsic calibration to temporal and spatial online integrated calibrations.

##  III METHODOLOGY

In this paper, a sensors pose estimation model is proposed to synchronize
sensors data from adjacent moments to the same moment. For the error term
equation of the pose estimation model, we adopted the correction based on the
environmental line feature to eliminate the influence of the error term on the
accuracy, due to the environmental line features are widely existed in indoor
and outdoor scenes and has a better correspondence between the camera and
LiDAR. Fig.2 illustrates the architecture of our method. Our proposed method
consists of four steps. First, the data between adjacent moments of the camera
and LiDAR are projected back to the same moment by means of pose estimation;
Second, the line features are extracted in a series of ways for both the
camera and LiDAR and are filtered and optimized; Then, a method for
eliminating dynamic point clouds is proposed, which only used the association
of adjacent two frames instead of detection boxed based on prior information.
Finally, the point cloud line features are projected on the pixel frames, and
the score is calculated and optimized during searching optimization periods.
Optimizing the score for calibrating the error of the previous stage of the
projection of the adjacent frames. The detailed steps are described in the
following sections.

###  III-A Pose estimation model

The temporal and spatial integrated calibration of camera and LiDAR lies in
the conversion matrix for temporal and spatial superposition. In this paper,
we assume that the camera acquires the image at moment $t$ and the LiDAR
acquires the point cloud at moment $t+e_{t}$ , the solution of the conversion
matrix between the adjacent moment is solved by the following method:
Estimating sensor motion within $e_{t}$ period and then calculating the
positional coordinate change between the vehicle body system and sensors
system onboard during $e_{t}$ period by means of coordinate system projection.
According to the pose estimation model:

$\begin{split}T&=T_{e}\times T_{v}\\\
T_{e}&=T_{e}^{\prime}+\epsilon_{e},T_{v}=T_{v}^{\prime}+\epsilon_{v}\end{split}$
(1)

where $T$ is the transform matrix from LiDAR at moment $t+e_{t}$ to camera at
moment $t$, $T_{e}$ is the measurement model of LiDAR-camera extrinsic
parameters, $T_{v}$ is the model of LiDAR-camera coordinate system motion
estimation as shown in Fig.3(a). Since there is an error between measurement
and ground truth, disassembly of $T_{e}$ and $T_{v}$ are shown, then Equation
1 is introduced in the next derivation:

$\begin{split}T=\left(T_{e}^{\prime}\cdot
T_{v}^{\prime}\right)+\left(T_{e}^{\prime}\epsilon_{v}+T_{v}^{\prime}\epsilon_{e}\right)+\epsilon_{e}\epsilon_{v}\end{split}$
(2)

where $T_{e}^{\prime}\cdot T_{v}^{\prime}$ is measured,
$\epsilon_{e}\epsilon_{v}$ is too small to be taken into account,
$T_{e}^{\prime}\epsilon_{v}+T_{v}^{\prime}\epsilon_{e}$ is the error.

Figure 3: A illustration of pose estimation model. (a) shows the pose
estimation model of the camera at time $t$ to the LiDAR at time $t+e_{t}$. (b)
shows the motion estimation via the inertial measurement unit navigation
conversion.

$T_{e}^{\prime}\epsilon_{v}$ represents the measurement error of the extrinsic
parameters, affected by initial values and prolonged vehicle driving.
$T_{v}^{\prime}\epsilon_{e}$ represents the error in the estimation of the
LiDAR-camera coordinate system motion. The motion of the inertial coordinate
system $T_{I_{-}I}$ during $e_{t}$ instant period is calculated by inertial
pre-integration as shown in Fig.3(b). Although the uniform linear motion
estimation model and the pre-integrated model is quite approximate to the
ground truth of the motion of the vehicle body coordinate system , the error
term has not to be eliminated since the motion of the LiDAR-camera coordinate
system is projected and estimated from the motion of the vehicle body
coordinate system, which reduces the accuracy of the projection between
adjacent moments $t$ and $t+e_{t}$. In order to correct the effect of the
error term online, line features are introduced in integrated calibration.

###  III-B Line features extraction

The line features are extracted to correct the error term in the previous
stage of optimization. As shown in Fig.LABEL:line_features_extraction, the
process of line feature extraction is carried out in both image and point
cloud. For image feature extraction, grayscale processing and canny
marginalization are adopted, then line feature extraction algorithm [20] is
applied. After that, the model of inverse distance transformation is
introduced to deal with the consistent correlation between LiDAR points and
image points, which leads to adopting a larger search step to prevent falling
into local optimum during searching optimization.

For LiDAR features extraction, the points are divided into different lines,
boundary line features are obtained by distance discontinuity.In order to
extract sufficient line features on low beam LiDAR, a local mapping method is
applied to combine three frames of point cloud into one, which can present
more points in one frame. Depending on the strength of the GPS signal and the
accuracy requirement of the autonomous driving scenario, we propose two local
mapping methods: the GPS-based approach and the Normal Distribution Transform
(NDT)-based approach[21]. The former method is less computationally intensive
with less accuracy, while the latter method is more computationally intensive
with more accuracy.

###  III-C Projection and feature filtering

For matching the line features extracted from the point cloud and the image,
pointcloud points set
$P^{L}=\left\\{p_{1}^{L},p_{2}^{L},p_{3}^{L},\cdots\cdots,p_{n}^{L}\right\\}$
and the image points set
$P^{C}=\left\\{p_{1}^{C},p_{2}^{C},p_{3}^{C},\cdots\cdots,p_{m}^{C}\right\\}$
are created, where $n$ is the number of point cloud points, $m$ is the number
of image points, $L$ is the LiDAR coordinate system, and $C$ is the camera
coordinate system. We define the rotation matrix $R_{L\rightarrow C}$ and the
translation vector $t_{L\rightarrow C}$ , which represent the coordinate
transformation process of the point cloud $P^{L}$ projected to $P^{C}$ in the
LiDAR coordinate system.

$P_{i}^{C}=R_{L\rightarrow C}\cdot P_{i}^{L}+t_{L\rightarrow C}$ (3)

Further filtering is done for outliers in line feature extraction. The above-
mentioned LiDAR point cloud has been converted into image form, so a
convolution kernel with an 8-pixel boundary is designed to filter out the
outliers on the image and get more organized line features from images. In
addition, for the points on the ground, are eliminated because their lateral
line features do not match well.

###  III-D Dynamic point clouds targets eliminating

Eliminating the dynamic point cloud targets is extremely essential to the
features matching. In the registration of point clouds at adjacent moments $t$
and $t+e_{t}$, although the stationary targets occupy the dominant position in
the observed environment, even a few dynamic targets will bring serious
interference to the accuracy of registration, especially in the case of
registration based on point cloud line feature rather than the registration
based on the original point cloud.

In order to remove dynamic targets, the current widely used method is target
detection based on deep learning. Origin point clouds and prior information as
the input of the detection network, which could identify possible dynamic
objects such as vehicles and pedestrians and return detection box. Then the
point clouds in the detection box are clustered and deleted to eliminate the
dynamic target. However, this method relies on prior information and an
additional detection network, which poses a real-time problem for online
integrated calibration.

In order to solve this problem, a lightweight method with excellent accuracy
has been proposed. The K-dimensional binary tree is used to quickly find the
nearest neighbors in the multidimensional space by using the dimensional
information of the segmentation. The nearest neighbor matching distance for
stationary objects is considered to be smaller, while the nearest neighbor
matching distance for dynamic objects is larger, so we set a threshold to
filter out dynamic targets. Because the point cloud at time $t$ is speculated
instead of ground truth of $t+e_{t}$, there should be an angular error.
According to the triangle similarity principle, the distant dynamic target
needs a larger threshold to be filtered out, so we set a factor to the fixed
threshold to make it a linear dynamic threshold. Dynamic targets are filtered
by linear dynamic thresholds and clustered for subsequent analysis.

Figure 4: An illustration of how the dynamic targets impact on line feature
projection. (a) is the result of projection without dynamic targets
eliminated, and (b) is the result of projection with dynamic targets
eliminated.  Algorithm 1 Optimization process

0: Image line features $I_{t}$ at frame $t$, Pointcloud horizontal
$F_{h}^{t+e_{t}}$ and vertical $F_{v}^{t+e_{t}}$ line features at frame
$t+e_{t}$, Initial extrinsic matrix $T_{e}$, Motion estimation matrix $T_{v}$,
last frame gray rate $gray\\_rate$

0: Calibrated extrinsic matrix and motion estimation;

1: Initialization: $score$, $max_{score}$ $\leftarrow$ 0\.

2: if $gray\\_rate\textgreater\gamma$ then

3: $step\\_size\\_larger=\alpha_{1}$

4: $step\\_size\\_smaller=\alpha_{2}$

5: else

6: $step\\_size\\_larger=\alpha_{3}$

7: $step\\_size\\_smaller=\alpha_{4}$

8: end if

9: for each $LiDAR\ point$ in $F_{h}^{t+e_{t}}$ do

10: $gray\\_value=weight*T_{e}*T_{v}*p_{t}$

11: $score+=gray\\_value$

12: end for

13: for each $LiDAR\ point$ in $F_{v}^{t+e_{t}}$ do

14: $gray\\_value=(1-weight)*T_{e}*T_{v}*p_{t}$

15: $score+=gray\\_value$

16: end for

17: if $score\textgreater max\\_score$ then

18: $max\\_score=score$

19: Update current_parameters

20: end if

21: $gray\\_rate=score/255/points\\_num$

22: return current_parameters

###  III-E Searching optimization

In the searching optimization process, an approach with both computation
accuracy and efficiency is introduced. In the previous stage, the LiDAR line
features have been extracted and projected onto the line features of the
image, and the proportion of LiDAR points projected onto the image gray area
has been calculated.

Figure 5: Different gray rate determines different step size strategies. (a)
shows a larger gray rate change and (b) shows a smaller gray rate change.

For computation accuracy, as is shown in Fig.5, four searching steps are
proposed according to different gray rates. For computation efficiency, a
searching method is applied to optimize the cost function. In [19], they
compare the current function score with adjacent 728 scores, if the searching
program finds parameters that have a higher score, it will stop the current
searching process and begin a new searching process at the position providing
a higher score. This searching process will stop when reaching the set
iteration count or finding the best score, thus being able to increase the
computation efficiency. We improved this process so that when the direction of
the fastest gradient descent is found, we follow that direction to get the
fastest optimization. As shown in Algorithm 1, this is the process of
searching optimization.

##  IV EXPERIMENT

In this section, the experiment setups and results of temporal and spatial
calibration are described in detail. We validate our proposed approach on
KITTI dataset[22], from which we used a high-resolution camera and a Velodyne
HDL-64E LiDAR. The LiDAR was scanned at 10 Hz and the data used in the
experiment is synced and rectified. The ground truth temporal and spatial
parameters can be obtained from calibration files.

###  IV-A Temporal calibration

Our method is compared to the approximate strategy of message filter in ROS.
Several evaluation metrics are designed in order to validate the accuracy,
which is based on Euclidean distance and ICP and NDT registration.

  * •

Euclidean distance-based evaluation metrics: The average Euclidean distance of
nearest point in the filter_cloud is calculated and compared.

  * •

NDT-based evaluation metrics: With the voxel downsampling grid size and other
parameters being the same, the final iteration number is calculated and
compared.

  * •

ICP-based evaluation metrics: with the maximum number of iterations, the
maximum distance and other parameters are the same, the score of Euclidean
distance is calculated and compared.

TABLE I: CALIBRATION RESULTS OF METRICS

Project | Dists_Evaluation_Score | NDT_Evaluation_Score | ICP_Evaluation_Score |
---|---|---|---|---
$infer\\_cloud\\_avg$ | 0.043 |  $stepsize1$ 4.091 |  $iterationNum1$ 0.042 |
$stepsize2$ 3.455 |  $iterationNum2$ 0.042 |
$infer\\_cloud\\_interpolation$ | 0.063 |  $stepsize1$ 19.10 |  $iterationNum1$ 0.054 |
$stepsize2$ 6.318 |  $iterationNum2$ 0.044 |
$filter\\_cloud\\_avg$ | 0.070 |  $stepsize1$ 33.82 |  $iterationNum1$ 0.059 |
$stepsize2$ 9.091 |  $iterationNum2$ 0.045 |

TABLE I shows the score three metrics: $infer\\_cloud\\_avg$ refers to the
average score of a point cloud treated by InferCloud at only 5Hz.
$infer\\_cloud\\_interpolation$ refers to the average score of a point cloud
interpolated upward at 10Hz. $filter\\_cloud\\_avg$ is the average score
processed using the ROS soft synchronization method. In NDT evaluation,
$stepsize1$ is 0.01m and $stepsize2$ is 0.05m. In ICP evaluation,
$iterationNum1$ is once while $iterationNum2$ is five. The whole results are
in Fig.6.

Figure 6: An illustration of three designed metrics performance. (a) shows the
Euclidean distance score of Euclidean distance-based evaluation metric; (b)
shows the final iteration number of NDT-based evaluation metrics; (c) shows
the Euclidean distance score of ICP-based evaluation metrics.

The data of Fig.6 tested based on the scenario where the point cloud frequency
is 5Hz, the image and inertial navigation frequency is 10Hz, and the point
cloud frame is 0.1s slower than the image. It reflects the advantages of our
method and the traditional ROS method in three evaluation metrics. In Fig.6,
the red lines ($filter\\_cloud$) represent the ROS soft synchronization score
curve, and the green and blue lines ($infer\\_cloud$) represent the score
curve after soft synchronization using our method. Meanwhile, in the NDT/ICP
evaluation metrics, we use different step sizes/different iterations to
reflect the universality of our method optimization. The green line represents
the smaller step length/fewer iterations, and the blue line represents the
larger step length/larger iterations.

Since the signal frequency processed by traditional ROS method is determined
by the data signal with the lowest frequency, it can be found that the
frequency of the $filter\\_cloud$ is 5Hz. If only the 5Hz point cloud is
processed by our method, and the solid green or blue line is compared with the
red line on the image, it can be found that our method is generally superior
to the ROS method in all aspects. In addition, our method can also interpolate
upward. The shallower line ($infer\\_cloud\\_interpolation$) below the
turquoise is the interpolation score curve. It can be seen that the extra
points processed by interpolation are still with the same level of
$filter\\_cloud$. Compared with the method of ROS, our approach improves 38.5%
in the Euclidean distance-based metrics and with excellent performance in NDT-
based evaluation metrics and ICP-based evaluation metrics.

###  IV-B Spatial calibration

We initially added a 2.0-degrees rotation disturbance on the X, Y, and Z axes
to the ground truth parameters. It must be clarified that whether the
2.0-degrees rotation disturbance is positive or negative is stochastic. During
the experiment, we compared the calibration error with the ground truth. In
addition, we tested the speed of correcting disturbance.

Figure 7: An illustration of instant correction of roll, pitch and yaw
disturbances. The red, green and blue lines in (a) represent the correction of
disturbance by our method after adding 2.0-degrees disturbance angle in roll,
pitch and yaw directions respectively. (b) shows the calibration result of our
method when the perturbation angle added 2.0-degrees in the roll direction.
(c) shows the calibration result of our method when the perturbation angle
added 2.0-degrees in the pitch direction. (d) shows the calibration result of
our method when the perturbation angle added 2.0-degrees in the yaw direction.

It can be seen from Fig.7 that our method can basically correct the extrinsic
parameter drift phenomenon caused by some reasons within 5 frames, and the
error generally shows a downward trend and can begin to converge at 1-3
frames. In Fig.7 (b),(c) and (d), the direction in which the deviation is
added is indicated by a darker color. It can be found that from the
perspective of the metric of angle correction, our method can also reduce the
error and basically complete convergence within 1-3 frames. But there is one
small problem that appeared on the direction of pitch angle correction,
speculated that for most of us choose picture perspective the middle section
of the line on the vertical direction features, and on the pitch direction, if
the deviation is also a small angle along with the vertical direction
deviation, such as the trunk of a tree after the middle period of line feature
extracting, even moving up and down will match the tree trunk well. The
overall impact of this oscillation problem is not very big, about $\pm 1.0$
degree. And it can be quickly corrected when there are obvious transverse line
features. The overall calibration results in more scenarios on the KITTI
dataset can be seen in Fig.8, which demonstrates that the proposed method is
applicable to different scenarios.

Figure 8: An illustration of spatial calibration of several scenarios. Green
points represent the projected LiDAR line points.

##  V CONCLUSION

The approach for online integrated temporal and spatial calibration for camera
and LiDAR is proposed. Based on the pose estimation model of the adjacent
moment, line features are selected to correct the error of the model. Hardware
interface and chessboard are not required in this integrated calibration
method. Besides, it’s more accurate than the soft simultaneous interpolation
method and upward interpolation between the point cloud between adjacent
moments is supported. It also demonstrates that the line features of point
clouds and images are robust features for matching and calibrating the error
of the pose estimation model. In addition, we show that the score of the
current integrated calibration results can be calculated and further exploited
to improve computational efficiency and accuracy. In future work, we plan to
introduce temporal and spatial online integrated calibration to multiple
cameras and multiple LiDARs.

## Acknowledgment

This paper and the research behind it would not have been possible without the
exceptional support of atatang Inc.(www.datatang.ai) that kindly provides the
High quality on demand professional services specialized in data collection
and annotation powering our model.

## References

  * [1] C. R. Qi, H. Su, K. Mo, and L. J. Guibas, “Pointnet: Deep learning on point sets for 3d classification and segmentation,” in Proceedings of the IEEE conference on computer vision and pattern recognition, pp. 652–660, 2017\.
  * [2] O. Ronneberger, P. Fischer, and T. Brox, “U-net: Convolutional networks for biomedical image segmentation,” in International Conference on Medical image computing and computer-assisted intervention, pp. 234–241, Springer, 2015\.
  * [3] T. Shan, B. Englot, D. Meyers, W. Wang, C. Ratti, and D. Rus, “Lio-sam: Tightly-coupled lidar inertial odometry via smoothing and mapping,” in 2020 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pp. 5135–5142, IEEE, 2020.
  * [4] J. Wu, Z. Xiong, et al., “A soft time synchronization framework for multi-sensors in autonomous localization and navigation,” in 2018 IEEE/ASME International Conference on Advanced Intelligent Mechatronics (AIM), pp. 694–699, IEEE, 2018.
  * [5] S. Kato, E. Takeuchi, Y. Ishiguro, Y. Ninomiya, K. Takeda, and T. Hamada, “An open approach to autonomous vehicles,” IEEE Micro, vol. 35, no. 6, pp. 60–68, 2015.
  * [6] A. Geiger, F. Moosmann, Ö. Car, and B. Schuster, “Automatic camera and range sensor calibration using a single shot,” in 2012 IEEE international conference on robotics and automation, pp. 3936–3943, IEEE, 2012\.
  * [7] J. Castorena, U. S. Kamilov, and P. T. Boufounos, “Autocalibration of lidar and optical cameras via edge alignment,” in 2016 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pp. 2862–2866, IEEE, 2016.
  * [8] A. English, P. Ross, D. Ball, B. Upcroft, and P. Corke, “Triggersync: A time synchronisation tool,” in 2015 IEEE International Conference on Robotics and Automation (ICRA), pp. 6220–6226, IEEE, 2015.
  * [9] P. H. Dana, “Global positioning system (gps) time dissemination for real-time applications,” Real-time systems, vol. 12, no. 1, pp. 9–40, 1997.
  * [10] M. Faizullin, A. Kornilova, and G. Ferrer, “Open-source lidar time synchronization system by mimicking gps-clock,” arXiv preprint arXiv:2107.02625, 2021.
  * [11] H. Sommer, R. Khanna, I. Gilitschenski, Z. Taylor, R. Siegwart, and J. Nieto, “A low-cost system for high-rate, high-accuracy temporal calibration for lidars and cameras,” in 2017 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pp. 2219–2226, IEEE, 2017.
  * [12] A. Harrison and P. Newman, “Ticsync: Knowing when things happened,” in 2011 IEEE International Conference on Robotics and Automation, pp. 356–363, IEEE, 2011.
  * [13] D. Scaramuzza, A. Harati, and R. Siegwart, “Extrinsic self calibration of a camera and a 3d laser range finder from natural scenes,” in 2007 IEEE/RSJ International Conference on Intelligent Robots and Systems, pp. 4164–4169, IEEE, 2007.
  * [14] M. Pereira, D. Silva, V. Santos, and P. Dias, “Self calibration of multiple lidars and cameras on autonomous vehicles,” Robotics and Autonomous Systems, vol. 83, pp. 326–337, 2016.
  * [15] Z. Taylor and J. Nieto, “Automatic calibration of lidar and camera images using normalized mutual information,” in Robotics and Automation (ICRA), 2013 IEEE International Conference on, Citeseer, 2013.
  * [16] R. Ishikawa, T. Oishi, and K. Ikeuchi, “Lidar and camera calibration using motions estimated by sensor fusion odometry,” in 2018 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pp. 7342–7349, IEEE, 2018.
  * [17] C.-L. Lee, Y.-H. Hsueh, C.-C. Wang, and W.-C. Lin, “Extrinsic and temporal calibration of automotive radar and 3d lidar,” in 2020 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pp. 9976–9983, IEEE, 2020.
  * [18] Z. Taylor and J. Nieto, “Motion-based calibration of multimodal sensor extrinsics and timing offset estimation,” IEEE Transactions on Robotics, vol. 32, no. 5, pp. 1215–1229, 2016.
  * [19] X. Zhang, S. Zhu, S. Guo, J. Li, and H. Liu, “Line-based automatic extrinsic calibration of lidar and camera,” in 2021 IEEE International Conference on Robotics and Automation (ICRA), pp. 9347–9353, IEEE, 2021.
  * [20] R. G. Von Gioi, J. Jakubowicz, J.-M. Morel, and G. Randall, “Lsd: a line segment detector,” Image Processing On Line, vol. 2, pp. 35–55, 2012.
  * [21] P. Biber and W. Straßer, “The normal distributions transform: A new approach to laser scan matching,” in Proceedings 2003 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS 2003)(Cat. No. 03CH37453), vol. 3, pp. 2743–2748, IEEE, 2003.
  * [22] A. Geiger, P. Lenz, C. Stiller, and R. Urtasun, “Vision meets robotics: The kitti dataset,” The International Journal of Robotics Research, vol. 32, no. 11, pp. 1231–1237, 2013.