Dictionary-Based Learning for 3D-Imaging with

Air-Coupled Ultrasonic Phased Arrays

Raphael Müller∗ , David Schenck∗ , Gianni Allevato† , Matthias Rutsch† , Jan Hinrichs† ,
Mario Kupnik† , and Marius Pesavento∗
∗ Communication Systems Group, Technische Universität Darmstadt, Germany
† Measurement and Sensor Technology Group, Technische Universität Darmstadt, Germany

Abstract—The Least Absolute Shrinkage and Selection Op-

erator (LASSO) was recently introduced for imaging of sparse
scenes in low sample size scenarios. In this paper, we investigate
the imaging capabilities of LASSO for ultrasonic phased arrays
operating in air. Using sparse regularization and a dictionary of
measured array point responses, 3D images can be generated with
only a few, single pulse measurements. A priori, the dictionary
of point responses is learned by calibrating the array. The
underlying low-rank structure can be adequately expressed by
fitting rank-one vectors. A proof of concept shows the imaging
potential of our proposed scheme with a real ultrasonic phased
array.
Index Terms—air-coupled ultrasound imaging, high resolution
imaging, array calibration, array characterization, dictionary
learning, best rank-one approximation, denoising compressed
sensing, sparse reconstruction, L1 regularization.

I. I NTRODUCTION Fig. 1. Uniform Rectangular Array (URA) consisting of 64 Murata MA404S

transducers. The array is controlled via custom electronics and features a
Ultrasonic imaging has a wide range of applications, e.g., 3D-printed waveguide for a grating-lobe free λ/2 acoustic aperture.
in navigation or healthcare [1], [2]. Compared to use cases in
fluids, airborne sensor solutions face a severe problem: The
II. M EASUREMENT S ETUP
reduced speed of sound in air and the requirement for large
coverage areas provoke large propagation delays which slow A. Air-Coupled Ultrasonic Phased Array
down imaging frame rates. This shortcoming motivates the Phased arrays are composed of a collection of spatially
use of sparse sampling concepts which provide high imaging distributed sensors with phase shifters that allow to electroni-
accuracy with only a few snapshot measurements. cally steer their receive and transmit beam. This enables both
Using a dictionary of array responses for different directions illumination of a distinct room section in the field of view
and ranges, the received array response to a transmitted pulse (FOV) of the sensor array and evaluation of the direction-of-
can be correlated to the entries of the dictionary in order to find arrival (DOA) for incoming point source signals. The amount
the corresponding angular and range locations present in the of sensors in the array and their inter-element distance define
received echo signal. However, for real hardware such array the array aperture, which gives a natural upper limit on the
responses are often difficult to characterize analytically. resolution capabilities of the system.
The focus of our research is therefore on the calibration of In our work, we consider an 8×8 uniform rectangular array
a given array in order to learn a dictionary of point responses (URA) consisting of ultrasonic transducers of type MA40S4S
based on sample measurements for said point responses. This manufactured by Murata (Nagaokakyō, Kyoto Prefecture,
approach is hardware-agnostic, making it applicable to various Japan). Key signatures of these transducers are an output
forms of ultrasonic imaging systems. The major contribution power level of more than 120 dB at a distance of 30 cm and
besides proposing an effective calibration procedure for air- their broad availability to both academic as well as commercial
coupled ultrasonic phased arrays is given by the successful customers. Their operating frequency of f0 = 40 kHz allows
imaging approach with real hardware based on a dictionary of us an imaging range of several meters as well as a sufficient
array responses. resolution of the scenery of interest.
This paper is structured as follows: In Section II, the
measurement setup is introduced. Section III describes the B. Customized Waveguide
LASSO imaging procedure followed by an explanation on how Achieving a grating-lobe free beam pattern is crucial to
to learn the dictionary in Section IV. Measurement results are prevent ghost targets during imaging. However, with a di-
given in Section V and Section VI concludes this paper. ameter of 9.9 mm the physical dimensions of the MA40S4S

978-1-7281-5448-0/20/$31.00 ©2020 IEEE

Authorized licensed use limited to: University of Canberra. Downloaded on May 23,2021 at 17:41:37 UTC from IEEE Xplore. Restrictions apply.
transducers do not allow to directly place them next to each position p = 1, . . . , P in the FOV. Having access to the angular
other at λ/2 ≈ 4.3 mm. In order to overcome this drawback, a response ap and frequency response bp , the array response dp
waveguide is used [3]. It creates an acoustic aperture with the is calculated via the Kronecker product
demanded λ/2-pitch. Furthermore, an essential benefit of this
setup is that the separation of the acoustic aperture from the dp = bp ⊗ ap . (3)
mechanical one offers the possibility to operate the array in Unfortunately, bp and especially ap are hard to describe
harmful environments, e.g., at elevated temperatures. A custom analytically, due to deviations in the signal reception of every
waveguide fabricated using standard 3D printing techniques is transducer and artifacts in the array geometry. Production
depicted in Fig. 1. tolerances and aging effects affect the phase and amplitude
Every waveguide channel adds the same path length from of received signals, thereby creating unequal sensor channels.
the acoustic aperture to the respective transducer. Whenever Customized arrays like the one at our hands which features an
the behavior of a single sensor or channel of the array is additional waveguide introduce further uncertainties into the
considered it should be noted that this includes both the point response behavior of the array. This motivates to empiri-
transducer and its respective waveguide channel. We therefore cally learn the array responses through sample measurements.
shift the focus from the mechanic to the acoustic aperture A calibration procedure to learn a dictionary from sample array
when referring to the phased array. The successful imaging responses is introduced next.
capabilities of the waveguide system have been proven earlier,
e.g., in [4], [5]. IV. L EARNING THE D ICTIONARY
III. I MAGING WITH S PARSE R EGULARIZATION A. Calibration
Assume that the scene of interest is overall sparse and only As shown in (3), the dictionary entry dp is constructed out
contains a few point targets. We have access to measurement of two - up to this point, unknown - array responses ap and
data received by M sensors, each one providing F frequency bp caused by a point source (point target) located at room
points covering the given signal bandwidth. Using a dictionary position p. The angular array response
D ∈ CM F ×P where P denotes the number of positions for
which a sample array response exists, imaging can be realized ap = a(ϑp , φp ) (4)
by applying the Least Absolute Shrinkage and Selection can be described as a function of azimuth ϑp and elevation
Operator (LASSO) on received signals. φp , whereas the frequency response
The LASSO procedure combines the classic Least Squares
(LS) problem with a sparsity promoting regularization term, bp = b(τp ) (5)
thereby introducing a trade-off between meeting the assump-
tion of a sparse scene and matching the receive data at can be described in dependence of the propagation delay τp
best to our dictionary. A minimization over the sum of a which is proportional to the range between the array and
data matching term under the Frobenius norm k · kF and a the target position p [5]. Note how bp is independent of
regularization term under the L1-norm k · k1 which sums up the transducer, thereby making the assumption of a common
the magnitudes of the entries in its argument expresses this frequency response for all sensors.
trade-off [6]: The active calibration procedure proposed in this paper
  involves an additional transducer placed at various room
min ||y − Dx||2F + µ||x||1 . (1) positions. It serves as an active echo target, transmitting pulses
x | {z } | {z } which are received by the phased array. At this point, the
Data matching Sparse regularization
imaging array is therefore only operated in its receive mode
The minimization problem in (1) aims at finding the solution to record the respective array responses. The corresponding
vector of possible target positions x ∈ CP ×1 out of P points measurement setup is depicted in Fig. 3.
in the FOV. Target locations are therefore obtained from the The total number of room positions P , sampled with this
support of the solution vector. In order to enforce sparsity on approach, not only defines the size of the dictionary D but
x, the regularizing parameter µ ∈ R weights the importance also the resolution of the image produced when solving the
of the L1-norm term. A large µ introduces a bias at the cost LASSO problem, as only discrete points are possible solutions
of ignoring the vectorized measurement snapshot y ∈ CM F ×1 to the optimization problem in (1). Hence, a fine grid of sample
obtained by illuminating the room. Consequently, a small value points is desirable for the calibration of D. Multiple snapshot
for the regularization parameter results in high confidence for measurements for every room point are considered in order to
the measurement data representing a sparse scene. further improve the precision of the measured array responses.
In general,
B. Low Rank Model
D = [d1 , d2 , . . . , dp , . . . , dP ] (2)
After capturing the received pulses emitted from the active
has the structure of a fat matrix (M F < P ) where each echo target at position p, the sampled time sequences of every
column dp represents the sample array response for one transducer are transformed to the frequency domain using the

Authorized licensed use limited to: University of Canberra. Downloaded on May 23,2021 at 17:41:37 UTC from IEEE Xplore. Restrictions apply.
Discrete Fourier Transform (DFT). Then, the DFT signals are
stored in a 3D-array

X p ∈ CM ×F ×L , (6)
3

where L denotes the number of snapshots used to determine

the point response vector for position p. 2
We assume that only the active transducer contributes to

ηp [%]
the echo signal received by each sensor element. Furthermore,
this active target is assumed to be a point source in the far 1

field of the array in an environment that is free of any multi-

path propagation. Consequently, if every transducer shows an 0
identical (yet potentially unknown) frequency response, then 60
the received signal in (6) is characterized by a low rank tensor 40
40.1
and can be decomposed into 40
m 20

X p = a(ϑp , φp ) ◦ b(τp ) ◦ hp (7) 39.9 f [kHz]

using the outer product ◦ of the array steering vector

Fig. 2. The relative reconstruction error ηp quantifies the model mismatch
a(ϑp , φp ) ∈ CM ×1 , the frequency response vector b(τp ) ∈ done when assuming the pre-processed X̃ X p to be of low rank. Most often,
CF ×1 and a complex pulse gain vector hp ∈ CL×1 [7]. some transducers m show a frequency response which deviates from the one
Since the dictionary entries in (3) depend on ap and bp assumed on basis of a low-rank model. Such deviations show explicitly at the
edges of the bandwidth under consideration.
only, the received data tensor X p in (6) is decomposed to
find the best rank-one approximations for the angular and
frequency responses related to position p. Therefore, a tensor D. Dictionary Evaluation
decomposition is applied which can be expressed in terms of
an LS minimization problem [6], [7] The accuracy of the tensor decomposition in (8) can be
evaluated using the relative reconstruction error
Xp − a
min kX b(ϑp , φp ) ◦ b(τ b p k2 ,
b p) ◦ h (8) L
F b p (`)|2
ap ,b
b p ,h
bp X p (m, f, `) − a
1 X |X̃ bp (m)bb p (f )h
ηp (m, f ) = (10)
b
L X p (m, f, `)|2
|X̃
`=1
where abp , b
b p and h b p are the estimates whose outer vector
product match the measured data X p best in terms of the for any signal channel m = 1, . . . , M and frequency point
Frobenius norm. Reconstructing the data tensor by the outer f = 1, . . . , F by considering all snapshots ` = 1, . . . L. The
product Xb p = abp ◦ b
bp ◦ h
b p also shows the benefit of denoising inspection of the relative error in (10) shows to which degree
the data as Xb p is free of noise components which cannot be the received, pre-processed data X̃ X p is of low tensor rank.
described using a low-rank structure. The relative error of the tensor decomposition for a specific
position within the FOV of the sensor is depicted in Fig. 2 for
C. Low-Rank Pre-Processing different frequencies and sensors.

Real ultrasonic transducers show a unique filter characteris- V. E XPERIMENTAL VALIDATION

tic that affects the spectrum of the received signal. However, In order to verify our imaging approach, we first calibrate
this violates our assumptions of a mutual frequency response our URA with an additional MA40S4S transducer acting as
equal for every transducer. Therefore, a pre-processing step is an active echo source. All measurements are conducted in an
proposed next that equalizes the data of every receive channel. anechoic chamber that features a goniometer with two radial
For every position p, the DFT signals stored in X p (m, f, `) and one linear axis to adjust ranges and directions with both,
are weighted by the respective absolute inverse value over high accuracy and repeatability, as shown in Fig. 3.
all snapshots and the absolute mean over all channels and We sample point responses for a FOV of ±20° in azimuth
snapshots, i.e., and elevation, respectively. In total, 1088 different angular
PL PM directions (ϑ, φ) at a distance of 2.1 m are measured with
X p (ν, f, ξ)|
ν=1 |X L = 4 snapshots for every room position p. This results in
X p (m, f, `) = X p (m, f, `) ξ=1
X̃ PL , (9)
M ξ=1 |X X p (m, f, ξ)| P = 1088 dictionary entries dp . Consequently, for each of
the P data tensors X p a rank-one approximation as shown in
for m = 1, . . . , M . (8) is performed.
The pre-processed X̃ X p is then subject to the proposed tensor Since our dictionary only contains point responses for a
decomposition in (8) and therefore replaces X p in the process single range, it is not possible to make any statement on the
of estimating a bp , b
b p and hbp. accuracy of the distance estimation of the echo target. Instead,

Authorized licensed use limited to: University of Canberra. Downloaded on May 23,2021 at 17:41:37 UTC from IEEE Xplore. Restrictions apply.
 1
20
0.9

0.8

10 0.7

0.6

φ [°]
0 0.5

0.4

−10 0.3

0.2

0.1
−20

Fig. 3. During calibration, a transducer (left) is mounted on a rail, firing 0

−20 −10 0 10 20
pulses at the phased array (right) which can be steered using two rotational ϑ [°]
axes. Together, rail and rotational axes allow to acquire point responses for
several ranges and directions.
Fig. 4. Generated 2D-image by solving the LASSO with a calibrated
dictionary for a FOV of ±20° based on a single snapshot measurement.
Plotting the solution |x| shows the correlation between the received echo
we now focus on a 2D-image for a proof of concept of the signal and corresponding dictionary entries. The scene only contains a single
proposed dictionary-based imaging approach. According to target, located at (ϑ, φ) = (10°, 10°) (marked in red).
our assumed model (3), the range and the direction information
of the response vector dp can be decoupled as described
by the frequency and directional array responses bp and ap , measurements allows to build the required dictionary from
respectively. Consequently, this allows separate imaging steps sampled receive data. This proves that the underlying rank-one
for the range and angle of targets relative to the phased array. model can indeed be applied to characterize a given ultrasonic
For the actual imaging trial, the additional transducer from phased array. Furthermore, our learning approach eliminates
the calibration step is mounted next to the URA to illuminate the need for analytic array responses which are tedious to de-
the room. The URA therefore still works in its receive mode scribe for real-world systems. Instead, our empirical approach
and the whole ultrasonic imaging system acts in a pitch- shows promising results for accurate imaging of sparse scenes.
catch setup. As a target, we place a single tetrahedral corner Future work focuses on validating our approach for a
reflector with an edge length of 14 cm at the sample distance full 3D-image, taking both direction and range into account.
2.1 m. Viewed from the array, the target is located at direction Additionally, resolving multiple targets is a subject of further
(ϑ, φ) = (10°, 10°). research.
In order to solve the LASSO, we make use of the Soft- R EFERENCES
Thresholding with Exact Line Search Algorithm (STELA)
[1] N. Gageik, P. Benz and S. Montenegro, “Obstacle Detection and Col-
[6]. This algorithm can be implemented on parallel hardware lision Avoidance for a UAV With Complementary Low-Cost Sensors,”
architectures to efficiently solve the optimization problem in in IEEE Access, vol. 3, pp. 599-609, May 2015.
(1). Hence, this contributes to our attempt on fast imaging [2] Y. Li et al., “High-speed integrated endoscopic photoacoustic and ultra-
sound imaging system,” IEEE Journal of Selected Topics in Quantum
cycles for air-coupled ultrasound. Electronics, vol. 25, no. 1, pp. 1-5, Jan.-Feb. 2019.
With the described setup, we are able to find a meaningful [3] A. Langen, Ein Verfahren zur Konstruktion anwendungsoptimierter Ul-
solution vector x with a single pulse measurement. In order traschallsensoren auf der Basis von Schallkanälen. IPA-IAO - Forschung
und Praxis. Heidelberg, Germany: Springer, 1993.
to achieve a sparse solution, the regularization parameter has [4] A. Jäger et al., “Air-coupled 40-kHz ultrasonic 2D-phased array based
been set to µ = 0.5kDT yk∞ where the term under the infinity on a 3D-printed waveguide structure,” IEEE International Ultrasonics
norm k · k∞ can be used as an initial value as a rule of thumb. Symposium (IUS), Washington, DC, Sept. 2017.
[5] G. Allevato et al., “3D imaging method for an air-coupled 40 kHz
Finally, Fig. 4 depicts the obtained 2D-image, where the ultrasound phased-array” 23rd International Congress on Acoustics
value of |x| provides an estimate on the position of the target. (ICA), Aachen, Germany, Sept. 2019.
The larger the value of |x| the more likely a target is positioned [6] Y. Yang, M. Pesavento, S. Chatzinotas and B. Ottersten, “Successive
convex approximation algorithms for sparse signal estimation with
at the corresponding grid point in the dictionary. nonconvex regularizations,” IEEE Journal of Selected Topics in Signal
Processing, vol. 12, no. 6, pp. 1286-1302, Dec. 2018.
VI. C ONCLUSION [7] G. Kushe, Y. Yang, C. Steffens and M. Pesavento, “A parallel sparse
regularization method for structured multilinear low-rank tensor decom-
Based on an experimentally measured dictionary of array position,” 27th European Signal Processing Conference (EUSIPCO), A
responses, we introduce a fast way to acquire images for air- Coruna, Spain, Sept. 2019.
coupled ultrasound phased arrays using sparse regularization.
Applying best rank-one approximation to array calibration

Authorized licensed use limited to: University of Canberra. Downloaded on May 23,2021 at 17:41:37 UTC from IEEE Xplore. Restrictions apply.