ARTIFICIAL LATERAL LINE SYSTEMS FOR FEEDBACK CONTROL OF
UNDERWATER ROBOTS
By
Ahmad Taha Abdulsadda

A DISSERTATION
Submitted to
Michigan State University
in partial fulﬁllment of the requirements
for the degree of
DOCTOR OF PHILOSOPHY
Electrical Engineering
2012

ABSTRACT
ARTIFICIAL LATERAL LINE SYSTEMS FOR FEEDBACK CONTROL OF
UNDERWATER ROBOTS
By
Ahmad Taha Abdulsadda
A lateral line system, consisting of arrays of ﬂow-sensing neuromasts, allows ﬁsh and
amphibians to probe their ambient environment and plays a vital role in their behaviors
spanning predator/prey detection, schooling, rheotaxis, courtship and communication. The
feats of biological lateral lines have inspired an increasing interest in engineering a similar
sensing module for underwater robots and vehicles. Often known as artiﬁcial lateral lines,
these sensors could potentially enable an underwater robot to detect and identify moving
or stationary objects, and exploit ambient ﬂow energy for eﬃcient locomotion. Despite the
advances made in this area, realizing a practical artiﬁcial lateral line still faces signiﬁcant
challenges in both signal processing and ﬂow sensor fabrication.
In this dissertation we describe our eﬀort in developing signal processing methods for
hydrodynamic object localization and tracking using an artiﬁcial lateral line (ALL). We
consider two types of objects, a vibrating sphere and a non-vibrating cylinder, both of which
are of interest in underwater applications. A vibrating sphere, known as a dipole source,
is widely used in the study of biological lateral lines and it emulates the rhythmic body or
ﬁn movement. A non-vibrating cylinder (with unknown cross-section shape), on the other
hand, represents a general moving or stationary object underwater with a 2D ﬂow proﬁle.
First, a novel bio-inspired artiﬁcial lateral line system is proposed for underwater robots
and vehicles by exploiting the inherent sensing capability of ionic polymermetal composites
(IPMCs). Analogous to its biological counterpart, the IPMC-based lateral line processes

the sensor signals through a neural network, and we demonstrate the performance of this
approach in the localization of a dipole source.
Second, with an assumption of potential ﬂows, we formulate nonlinear estimation problems for the localization and tracking of a dipole source based on analytical ﬂow models, and
propose and compare several algorithms for solving the problem. For the case of a moving
cylinder, we use conformal mapping to represent a general cross-section proﬁle, and explore
the use of Kalman-ﬁltering-type techniques in the tracking and size/shape estimation of the
object.
We have conducted extensive experiments to validate the developed algorithms with an
artiﬁcial lateral line prototype made of millimeter-scale IPMC sensors, with sensor-to-sensor
separation of 2 cm, which is determined through an optimization process based on the
Cramer-Rao bound (CRB) analysis.
Finally, we experimentally explore the use of IPMC sensors for estimating the hydrodynamic parameters involved in a K´rm´n vortex street that is created by a stationary cylinder
a a
in a ﬂow. We validate that the vortex shedding frequency, which can be extracted from the
sensor signal, shows clear correlation with the ﬂow speed and the obstacle size.

Copyright by
AHMAD TAHA ABDULSADDA
2012

DEDICATION

To Mom, Dad, and My Family

v

ACKNOWLEDGMENTS

I am grateful to my advisor, Dr. Xiaobo Tan, for his thoughtful guidance during the course
of my PhD study. He inspired me with his expertise and guided me with his patience in the
novel ﬁeld of artiﬁcial lateral lines for feedback control of underwater robots.
I would like to thank my committee members, Dr. Hassan Khalil, Dr. Fathi Salem, and
Dr. Guoming (George) Zhu for providing helpful suggestions and comments. Thanks also
go to my colleagues for their help and discussion, in particular, John Thon, Hong Lie, and
Fetain Zhang.
I am grateful for the ﬁnancial support of my study from Ministry of Higher Education
in Iraq and my research from the Oﬃce of Naval Research (ONR: Grant N000140810640,
N000141210149) and National Science Foundation (NSF: ECCS 0547131, CCF 0820220, CNS
0751155, DBI-0939454, IIS 0916720).
I would like to thank all members in Iraqi Cultural Oﬃce in Washington DC. Finally,
I would like to thank my family for their patience, understanding and support during this
study.

vi

TABLE OF CONTENTS

List of Tables . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

x

. . . . . . . . . . . . . . . . . .

xii

Chapter 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.1 Research Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2 Research Objective . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3 Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3.1 An Artiﬁcial Lateral Line System Using IPMC Sensor Arrays . . . .
1.3.2 Localization of a Fixed Dipole Using Neural Network . . . . . . . . .
1.3.3 Nonlinear Estimation-based Dipole Source Localization . . . . . . . .
1.3.4 Localization of a Moving Dipole Source Underwater . . . . . . . . . .
1.3.5 Underwater Tracking and Size-Estimation of a Moving Object . . . .
1.3.6 Localization and Identiﬁcation of an Stationary Obstacle Underwater
1.4 Organization of This Dissertation . . . . . . . . . . . . . . . . . . . . . . . .

1
1
5
6
6
8
8
9
10
11
12

Chapter 2 Dipole Source Localization Using Neural Networks .
2.1 Hydrodynamic Model for an IPMC Sensor Under Oscillatory Flow
2.2 Source Localization Using Neural Network Processing . . . . . . .
2.3 Optimum Design of the Neural Network Using Genetic Algorithm
2.4 Experimental Setup and Sensor Characterization . . . . . . . . .
2.5 Experimental Results on Diploe Source Localization . . . . . . . .
2.6 Chapter Summary . . . . . . . . . . . . . . . . . . . . . . . . . .

List of Figures

. . . . . . . . . . . . . .

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

14
15
19
22
23
29
32

Chapter 3 Nonlinear Estimation-based Dipole Source Localization
3.1 Problem Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2 Estimation Algorithms . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2.1 Gauss-Newton Algorithm . . . . . . . . . . . . . . . . . . . . .
3.2.2 Newton-Raphson Algorithm . . . . . . . . . . . . . . . . . . .
3.2.3 Beamforming Algorithm . . . . . . . . . . . . . . . . . . . . .
3.2.4 Cramer-Rao Bound Analysis . . . . . . . . . . . . . . . . . . .
3.3 Simulation Results . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.1 CRB-based Optimal Design of Lateral Line Sensor . . . . . . .

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

37
38
42
43
44
45
46
49
49

vii

.
.
.
.
.
.
.

3.3.2

3.4

3.5

Localization of Dipole Source with Unknown Vibration Amplitude and
Orientation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.3 Further CRB Analysis on Performance . . . . . . . . . . . . . . . . .
Experimental Setup and Flow Model Veriﬁcation . . . . . . . . . . . . . . .
3.4.1 Experimental Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.4.2 Flow Model Validation . . . . . . . . . . . . . . . . . . . . . . . . . .
3.4.3 Calibration of IPMC Sensors . . . . . . . . . . . . . . . . . . . . . . .
3.4.4 Experimental Results on Dipole Source Localization . . . . . . . . . .
Chapter Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Chapter 4 Localization of a Moving Dipole Source . . . . . . . . .
4.1 Localization of a Traveling Dipole Source: Problem Formulation . .
4.2 The Sliding Discrete Fourier Transform (SDFT) Algorithm . . . . .
4.3 The Gauss-Newton Method for Solving the Estimation Problem . .
4.4 Simulation Results . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.5 Experimental Setup and Results . . . . . . . . . . . . . . . . . . . .
4.5.1 Experimental Setup . . . . . . . . . . . . . . . . . . . . . . .
4.5.2 Experimental Results on Moving Dipole Source Localization
4.6 Chapter Summary . . . . . . . . . . . . . . . . . . . . . . . . . . .
Chapter 5 Underwater Tracking and Size-Estimation
5.1 Flow Model and Problem Formulation . . . . . . . .
5.2 Nonlinear Filtering . . . . . . . . . . . . . . . . . . .
5.2.1 Extended Kalman Filter . . . . . . . . . . . .
5.2.2 Unscented Kalman Filter . . . . . . . . . . . .
5.2.3 Particle Filter . . . . . . . . . . . . . . . . . .
5.3 Simulation Results . . . . . . . . . . . . . . . . . . .
5.3.1 Recursive Information Matrix . . . . . . . . .
5.3.2 Simulation Results on Nonlinear Filtering . .
5.4 Experimental Results . . . . . . . . . . . . . . . . . .
5.5 Chapter Summary . . . . . . . . . . . . . . . . . . .

of a
. . .
. . .
. . .
. . .
. . .
. . .
. . .
. . .
. . .
. . .

Chapter 6 Exploration of IPMC Sensors in Vortex Street
6.1 Stationary Vortex and K´rm´n Vortex Street . . . . . . .
a a
6.2 Localization of an Obstacle: Problem Formulation . . . . .
6.3 Experimental Setup and Results . . . . . . . . . . . . . . .
6.3.1 Swim Tunnel Flow Calibration . . . . . . . . . . . .
6.3.2 Diﬀerentiation of Steady and Unsteady Flow . . . .
6.3.3 Diﬀerentiation of Obstacle Size . . . . . . . . . . .
6.3.4 Diﬀerentiation of Flow Velocity . . . . . . . . . . .
6.3.5 Two IPMC Sensors . . . . . . . . . . . . . . . . . .
6.4 Chapter Summary . . . . . . . . . . . . . . . . . . . . . .
Chapter 7

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

72
73
76
78
79
83
83
86
87

Moving Object
. . . . . . . . . .
. . . . . . . . . .
. . . . . . . . . .
. . . . . . . . . .
. . . . . . . . . .
. . . . . . . . . .
. . . . . . . . . .
. . . . . . . . . .
. . . . . . . . . .
. . . . . . . . . .

96
97
101
102
103
105
109
109
111
115
118

Sensing
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .

127
128
133
137
137
139
142
143
143
144

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

51
53
58
58
59
61
64
70

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152
viii

Chapter 8

Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155

Appendices . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . 158

Appendix A Artiﬁcial Lateral Line Design: CRB

. . . . . . . . . . . . . . . . 159

Appendix B Appendix for Chapter5 . . . . . . . . . . . . . . . . . . . . . . . . . 161
B.1 Derivation of CRB . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 161
B.2 More Results for One-stage Extended Kalman Filtering Scheme . . . . . . . 164
Bibliography . . . . . . . . . . . . . . .

ix

. . . . . . . . . . . . . . . . . . . 169

LIST OF TABLES

Table 3.1

Table 3.2

Experimental results: maximum and average estimation errors under diﬀerent algorithms. “BF (0.02)” stands for beamforming with
0.02×0.02 cm2 grid, and similar interpretations apply for “BF (0.05)”
and “BF (0.1)” . Note that vibration amplitude and orientation are
assumed known for the BF algorithms. . . . . . . . . . . . . . . . .

68

Average computational time per localization point for diﬀerent algorithms. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

70

Table 4.1

Simulation results: Estimated velocity (vm ), vibrating amplitude (A). 85

Table 4.2

Experimental results: Estimated velocity (vm ), vibrating amplitude
(A). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

87

Table 5.1

Simulation results: Estimated velocity, size parameter, and shape parameter under diﬀerent process noise variances. Two-stage extended
Kalman ﬁltering is used. . . . . . . . . . . . . . . . . . . . . . . . . 114

Table 5.2

Simulation results: Estimated velocity, size parameter, and shape
parameter under diﬀerent measurement noise variances. Estimation
is done with two-stage extended Kalman ﬁltering. . . . . . . . . . . 114

Table 5.3

Simulation results: Estimated velocity, size parameter, and shape
parameter under three nonlinear ﬁltering algorithms. . . . . . . . . . 115

Table 5.4

Cylindrical object experimental results: Estimated velocity and size
parameter. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117

Table 5.5

Ellipsoidal object experimental results: Estimated velocity, size parameter, and shape parameter. . . . . . . . . . . . . . . . . . . . . . 117

Table 6.1

Simulation results: Setup parameters. . . . . . . . . . . . . . . . . . 132
x

Table B.1

One-stage scheme: Estimated velocity, size parameter, and shape
parameter under diﬀerent process noise variances. . . . . . . . . . . 164

Table B.2

One-stage scheme: Estimated velocity, size parameter, and shape
parameter under diﬀerent measurement noise variances. . . . . . . . 167

xi

LIST OF FIGURES

Figure 1.1

(a) Illustration of the structure of a neuromast (image credit: C.
H. Mallery); (b) distribution of superﬁcial (small dots) and canal
neuromasts (dots within shaded areas) on the Lake Michigan mottled
sculpin [20]. For interpretation of the references to color in this and
all other ﬁgures, the reader is referred to the electronic version of this
dissertation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3

Illustration ot the potential use of an artiﬁcial lateral line system on
a robotic ﬁsh. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6

Figure 1.3

Illustration of IPMC sensing principle. . . . . . . . . . . . . . . . . .

7

Figure 1.4

Beneﬁt of exploiting vortex streets (adapted from [57]): (a) Exploiting
ﬂow due to vortex streets of leading ﬁsh; (b) harvesting energy from
oncoming vortices. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

12

Figure 2.1

Illustration of the beam dynamics of an IPMC sensor in ﬂow. . . . .

16

Figure 2.2

Schematic of the dipole source localization problem. . . . . . . . . .

19

Figure 2.3

Schematic of the MLP neural network for signal processing of the
IPMC lateral line. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

20

Figure 2.4

A prototype of IPMC-based lateral line system. . . . . . . . . . . . .

23

Figure 2.5

The schematic of circuit for measuring the short-circuit current generated by the IPMC sensor. . . . . . . . . . . . . . . . . . . . . . .

24

Experimental setup: (a) Experimental tank equipped with a DPIV
system; (b) the mini-shaker-based dipole source. . . . . . . . . . . .

25

Figure 1.2

Figure 2.6

xii

Figure 2.7

Figure 2.8

Figure 2.9

Figure 2.10

Figure 2.11

Figure 2.12

Figure 2.13

Figure 2.14

Figure 2.15

Figure 3.1

Figure 3.2

Figure 3.3

Figure 3.4

Flow ﬁeld around one IPMC sensor under a dipole stimulation: (a)
Dipole-sensor separation 10 cm, with maximum ﬂow amplitude of 9.3
mm/s in the ﬁeld; (b) dipole-sensor separation 45 cm, with maximum
ﬂow amplitude of 4.1 mm/s in the ﬁeld. . . . . . . . . . . . . . . . .

26

Typical IPMC sensor signals: (a) Signal from sensor #2, which is
close to the dipole source; (b) signal from sensor #4. . . . . . . . . .

27

Sensor characterization, measured sensor signal magnitude as a function of source location, The dipole vibration axis is: (a) Parallel to
the IPMC beam plane; (b) perpendicular to the IPMC beam plane.

28

Evolution of the ﬁtness function (the mean square error J) over generations, in the genetic algorithm-based training phase. This example
involves all six sensors and 180 training points. . . . . . . . . . . .

30

Localization of the dipole source using all six sensors of the lateral
line: (a) With 90 training points; (b) with 50 training points. . . . .

31

Localization error at 16 validation points when all six sensors are
used: (a) With 90 training points; (b) with 50 training points. . . .

32

Localization of the dipole source using four sensors of the lateral line:
(a) With 90 training points; (b) with 50 training points. . . . . . . .

34

Localization of the dipole source using two sensors of the lateral line:
(a) With 90 training points; (b) with 50 training points. . . . . . . .

35

Error in localization of the dipole source using the lateral line: (a)
Maximum error along the track; (b) average error along the track. .

36

Illustration of the problem setup. The dipole source is located at
(xs , ys ). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

40

CRB-based cumulative uncertainty in localization as a function of
intra-sensor spacing. . . . . . . . . . . . . . . . . . . . . . . . . . . .

51

Illustration of the simulation setup, where 19 points as indexed in the
ﬁgure will be localized. . . . . . . . . . . . . . . . . . . . . . . . . .

53

Simulation results: (a) Estimation of source locations; (b) localization
errors at the 19 points. . . . . . . . . . . . . . . . . . . . . . . . . .

54

xiii

Figure 3.5

Simulation results: Estimation of the source vibration amplitudes. .

55

Figure 3.6

Simulation results: Estimation of the source vibration orientations. .

56

Figure 3.7

Simulation results: Comparison of square-root of empirical estimation
variances (a) σx and (b) σy with the corresponding CRBs. . . . . .

57

Figure 3.8

Experimental setup: Schematic of the experimental system. . . . . .

60

Figure 3.9

Experimental results: Comparison of DPIV-measured ﬂow velocity
with the model prediction. . . . . . . . . . . . . . . . . . . . . . . .

61

Figure 3.10

Illustration of IPMC sensor calibration.

. . . . . . . . . . . . . . .

63

Figure 3.11

Experimental results: (a) Estimation of source locations; (b) localization errors at the 19 points. . . . . . . . . . . . . . . . . . . . . .

65

Experimental results: (a) Estimation of the source vibration amplitudes; (b) estimation of the source vibration orientations. . . . . . .

66

Figure 3.12

Figure 3.13

Figure 3.14

Experimental results: Constructed energy-like map when the dipole
is located at point #3 (0.1 × 0.1 cm2 grid). . . . . . . . . . . . . . .

67

Experimental results: Localization of the dipole source using the
beamforming algorithm: (a) Estimation of the source locations; (b)
localization errors at the 19 points. . . . . . . . . . . . . . . . . . .

69

Illustration of the problem setup for the localization of a moving
dipole source. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

74

Figure 4.2

Illustration of the simulation setup. . . . . . . . . . . . . . . . . . .

80

Figure 4.3

Simulation results: Time trajectories of the ﬂow velocities at the
sensor locations. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

81

Simulation results: Magnitude trajectories of the DC component at
the sensor sites. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

82

Simulation results: Magnitude trajectories of the AC component at
the sensor sites. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

83

Figure 4.1

Figure 4.4

Figure 4.5

xiv

Figure 4.6

Simulation results on the estimation of the source location: (a) Localization; (b) tracking error. . . . . . . . . . . . . . . . . . . . . . .

84

Figure 4.7

Simulation results on the estimation of the source location, vibration
amplitude changing with the time: (a) Localization; (b) tracking error. 88

Figure 4.8

Simulation results of the source vibration amplitude. . . . . . . . . .

89

Figure 4.9

An experimental prototype of IPMC-based lateral line system. . . .

89

Figure 4.10

Experimental setup of the moving dipole source: Schematic. . . . . .

90

Figure 4.11

Experimental setup of the moving dipole source: Picture of the experimental setup. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

91

Figure 4.12

Calibration of an IPMC sensor.

. . . . . . . . . . . . . . . . . . . .

91

Figure 4.13

Experimental results: Time trajectories of the ﬂow velocities at the
sensor locations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

92

Experimental results: Magnitude trajectories of the DC component
at the sensor sites. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

93

Experimental results: Magnitude trajectories of the AC component
at the sensor sites. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

94

Experimental results on the estimation of the source location: (a)
Localization; (b) tracking error. . . . . . . . . . . . . . . . . . . . .

95

Illustration of the Laurent series expansion : (a) Circle ζ(z) = z;
(b) the shape ﬁngerprint ζ(z) = z + 1/(2z); (c) the shape ﬁngerprint
ζ(z) = z + 1/(2z) + 1/(6z 2 ); (d) the shape ﬁngerprint ζ(z) = z +
1/(2z) + 1/(6z 2 ) + 1/(12z 3 )(adapted from [11]). . . . . . . . . . . .

99

Figure 4.14

Figure 4.15

Figure 4.16

Figure 5.1

Figure 5.2

Illustration of the problem setup for the tracking and estimation of a
moving object. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100

Figure 5.3

Illustration of the two-stage ﬁltering scheme. . . . . . . . . . . . . . 108

Figure 5.4

√

CRB for: (a) Moving object xs coordinate; (b) moving object ys
coordinate; (c) size parameter R; (d) shape parameter λ1 . . . . . . . 110
xv

Figure 5.5

Simulation setup for tracking of a moving cylinder. . . . . . . . . . . 111

Figure 5.6

Simulated time trajectories of sensor measurements. . . . . . . . . . 112

Figure 5.7

Simulation results on tracking the ellipsoidal cylinder using two-stage
extended Kalman ﬁltering: (a) The estimation of the moving object
location for diﬀerent variance values for the process noises; (b) localization error at each x-coordinate of the object. The measurement
noise variance is held at 0.01 cm2 /s2 . . . . . . . . . . . . . . . . . . 113

Figure 5.8

Simulation results on tracking the ellipsoidal cylinder using the twostage extended Kalman ﬁltering: (a) The estimation of the moving object location for diﬀerent variance values for the measurement
noises; (b) localization error at each x-coordinate of the object. The
process noise variance is held at 0.3 (in its respective unit for each
state component). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119

Figure 5.9

Simulation results on tracking the ellipsoidal cylinder for three nonlinear ﬁlters: The process noise variance is held at 0.3 (in its respective
unit for each state component) and the measurement noise variance
is held at 0.01 cm2 /s2 . . . . . . . . . . . . . . . . . . . . . . . . . . 120

Figure 5.10

Experimental setup: Schematic. . . . . . . . . . . . . . . . . . . . . 121

Figure 5.11

Experimental setup: (a) Picture of the experimental setup with a
circular cylinder; (b) picture of the experimental setup with an ellipsoidal cylinder. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122

Figure 5.12

Comparison of DPIV-measured ﬂow velocity with the model prediction.123

Figure 5.13

Measured ﬂow velocities by six IPMC sensors: (a) DC motor driving
the cylindrical object, with digital low pass ﬁlter (fc = 42) Hz; (b) AC
motor driving ellipsoidal object, with digital low pass ﬁlter (fc = 4) Hz.124

Figure 5.14

Experimental results on tracking a cylinder with circular proﬁle with
two-stage extended Kalman ﬁltering: (a) The estimation of the object
location; (b) localization error at each x−coordinate of the object. . 125

Figure 5.15

Ellipsoidal object experimental results on tracking an ellipse: (a)
The estimation of the object location; (b) localization error at each
x−coordinate of the object. . . . . . . . . . . . . . . . . . . . . . . . 126
xvi

Figure 6.1

Illustration of the analytical model for cylinder-shed vorticity. . . . . 129

Figure 6.2

Visualization of the ﬂow around a cylinder. . . . . . . . . . . . . . . 130

Figure 6.3

Illustration of the analytical model for cylinder-shed K´rm´n vortex
a a
street. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130

Figure 6.4

Visualization of a K´rm´n vortex street ﬂow ﬁled around: (a) Circular
a a
cylinder; (b) ellipsoidal cylinder. . . . . . . . . . . . . . . . . . . . . 133

Figure 6.5

Illustration of the problem setup. The obstacle is located at (xo , yo ). 134

Figure 6.6

Simulation results: 3D illustration of ﬂow velocity behind the obstacles: (a) Along the x−direction; (b) along the y−direction. In each
subﬁgure, the left plot shows the velocity for the circular cylinder
while the right plot shows the velocity for the ellipsoidal cylinder. . . 136

Figure 6.7

Loligo Systems swim tunnel: (a) Side view; (b) top view. . . . . . . 138

Figure 6.8

Rotor-X ﬂow sensor for swim tunnel calibration: (a) Side view; (b)
top view; (c) calibration curve . . . . . . . . . . . . . . . . . . . . . 139

Figure 6.9

IPMC sensor in a uniform ﬂow in the absence of the obstacle: (a)
Sensor at to ; (b) sensor at t1 , (b) sensor at t2 . . . . . . . . . . . . . 141

Figure 6.10

IPMC sensor right behind obstacle (predominantly negative ﬂow):
(a) Illustration; (b) typical sensor conﬁguration. . . . . . . . . . . . 141

Figure 6.11

IPMC sensor behind obstacle but shifted from center (predominantly
positive ﬂow): (a) Illustration; (b) typical sensor conﬁguration. . . . 142

Figure 6.12

Typical IPMC sensor signals and frequency spectrum of the response
under diﬀerent ﬂow conditions with nominal ﬂow velocity of 10.6 cm/s:
(a) Running water (free ﬂow); (b) predominantly negative ﬂow (suction region), frequency = 0.61 Hz; (c) predominantly positive ﬂow
(vortex region), frequency = 0.67 Hz. . . . . . . . . . . . . . . . . . 145

Figure 6.13

IPMC sensor was attached behind of the: (a) 6 cm diameter cylindrical obstacle; (b) 4 cm diameter cylindrical obstacle. . . . . . . . . 146
xvii

Figure 6.14

Typical IPMC sensor signals for diﬀerent-size cylindrical obstacles
with diameter: (a) 6 cm, theoretical / measured fsh : 0.56/0.6 Hz;
(b) 4 cm, theoretical / measured fsh : 0.73/0.79 Hz. . . . . . . . . . 147

Figure 6.15

Typical IPMC sensor signals under 6 cm upstream cylinder, the nominal ﬂow velocity of: (a) 7.5 cm/s, Theoretical / Measured fsh :0.4/0.61
Hz;(b) 10.6 cm/s, Theoretical / Measured fsh : 0.56/0.67 Hz. . . . . 148

Figure 6.16

Typical IPMC sensor signals to the human ﬁngers that moved alongside the two IPMC sensors with velocity of about 25 cm/s, direction
of the movement: (a) Toward positive x-axis; (b) toward negative
x-axis. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 149

Figure 6.17

Two IPMC sensors were attached behind a 6 cm cylinder: (a) Setup
at time to ; (b) sensor #1 sucked toward the obstacle, while sensor#2
repelled from the obstacle at time t1 . . . . . . . . . . . . . . . . . . 150

Figure 6.18

Signals from the two IPMC sensors and their frequency spectra: (a)
Sensor #1, dominant frequency = 0.55 Hz; (b) sensor #2, dominant
frequency = 0.73 Hz. . . . . . . . . . . . . . . . . . . . . . . . . . . 151

Figure 8.1

Illustration of the n-stages nonlinear ﬁltering scheme. . . . . . . . . 156

Figure B.1

Simulation results on tracking the ellipsoidal cylinder for one-stage
scheme: (a) The estimation of the moving object location for diﬀerent
variance values for the process noises; (b) localization error at each
x-coordinate of the object. The measurement noise variance is held
at 0.01 cm2 /s2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 165

Figure B.2

Simulation results on tracking the ellipsoidal cylinder for one-stage
scheme: (a) The estimation of the moving object location for diﬀerent
variance values for the measurement noises; (b) localization error at
each x-coordinate of the object. The process noise variance is held at
0.3 (in its respective unit for each state component). . . . . . . . . . 166

xviii

Chapter 1
Introduction

1.1

Research Background

Most ﬁsh and aquatic amphibians use the lateral line system as an important sensory organ
to probe their environment [20, 21]. A lateral line consists of arrays of hair cell sensors,
known as neuromasts. Each neuromast contains bundles of sensory hairs, encapsulated in a
gelatinous structure called cupula, as illustrated in Fig. 1.1(a). An impinging ﬂow deﬂects
the cupula, and thus the hairs inside, eliciting ﬁring of the hair cell neurons. Neuromasts can
be divided into two types, superﬁcial neuromasts, which are distributed on the skin surface,
and canal neuromasts, which are recessed in the scales or in bony canals underneath the skin
(Fig. 1.1(b)). With the same basic structure, the two types of neuromasts show distinct sensing characteristics [52]. The lateral line system allows an aquatic animal to identify near-ﬁeld
objects of interest and perform hydrodynamic imaging of the environment, typically within
one to two Body Lengths (BLs) of the animal. Consequently, the lateral line is involved in
various behaviors of aquatic animals, such as prey/predator detection [70], schooling [69],
1

rheotaxis [36], courtship and communication [21]. In addition to the qualitative roles the
lateral line plays in behaviors, there have been studies on how probed information is encoded
and decoded in the nervous systems [26]. For example, Fishes frequently exploit unsteady
forcing of wake ﬂows to extract energy for eﬃcient locomotion and maneuvering. The researchers in [56, 55] observed that, when trout perform K´rm´n gaiting to hold station in
a a
the wake of vortices shed from a cylindrical obstacle, their muscle activity is signiﬁcantly
reduced comparing to swimming in steady ﬂow. The role of the lateral line organ in such
behavior is further described in [93].
An engineering equivalent of a biological lateral line is of great interest to the navigation
and control of underwater robots and vehicles. In particular, an artiﬁcial lateral line will represent a new, noiseless sensing modality for underwater applications that is complementary
to traditional sensors such as vision and sonar [20, 59]. A number of micromachined ﬂow
sensors, inspired by neuromasts in ﬁsh lateral lines [74, 31, 96, 60] or by wind receptor hairs
in insects [67, 29], have been reported in the last decade. These sensors typically have an
out-of-plane beam structure that bends or deforms during the interaction with a ﬂow. The
bending or deformation, which carries information about the ﬂow, is captured via resistance
change [67, 92, 74, 31, 60, 97] or capacitance change [78, 29]. Hair cell-inspired sensors have
also been developed at slightly larger scales, exploiting optical transduction [50] or novel
sensing materials such as ionic polymer-metal composites (IPMCs) [1, 2] and gel-supported
lipid bilayers [80].
Extensive behavioral and neurophysiological studies have been conducted by biologists to
understand how hydrodynamic stimuli are encoded and processed in the lateral line system
[26, 24, 94, 9, 22, 23], which could potentially provide insight into information processing
2

(a)

(b)
Figure 1.1: (a) Illustration of the structure of a neuromast (image credit: C. H. Mallery);
(b) distribution of superﬁcial (small dots) and canal neuromasts (dots within shaded areas)
on the Lake Michigan mottled sculpin [20]. For interpretation of the references to color in
this and all other ﬁgures, the reader is referred to the electronic version of this dissertation.

for artiﬁcial lateral lines. Theoretical work on ﬂow modeling in the context of lateral lines
is often instrumental in both explaining hydrodynamic imaging by ﬁsh and guiding the
information extraction in artiﬁcial lateral lines. For example, Hassan modeled the ﬂow ﬁeld
around a ﬁsh when the latter passes by or approaches an obstacle [40], moves in open water
or approaches a plane surface [41], glides alongside or above a plane surface [42], and is near
an oscillating sphere (also known as dipole) [43]. Ren and Mohseni derived the pressure
3

ﬁeld around a ﬂat body as well as the ﬂow inside a trunk canal in the presence of vortex
street and explored its implications in vortex sensing [75]. Sichert et al. explored the use
of hydrodynamic multipole expansion to explain how aquatic animals distinguish shapes of
moving objects [83], and in another work Sichert and van Hemmen investigated the shape
inﬂuence of a submerged moving object on the lateral line perception using the conformal
mapping theory [84]. Bouﬀanais et al. performed analytical and numerical studies on the
pressure ﬁeld around a cylindrical object in a uniform ﬂow and examined the estimation of
location, size, and cross-section shape of the object based on the measured pressure ﬁeld
[11]. There has also been experimental work in exploring the strategies for vortex sensing
[50, 4, 91] and object tracking and recognition [33].
Barring some of the aforementioned research in vortex sensing and object imaging, much
of the work in information processing for biological and artiﬁcial lateral lines has been focused on dipole source localization. The dipole source emulates the rhythmic movement of
ﬁns and body appendages, and has been widely used as a stimulus (playing the role of predator, prey, or conspeciﬁc) in the study of biological lateral lines [26, 22, 23, 39]. Dipole source
localization has also become a benchmark problem in the development of artiﬁcial lateral
lines, for demonstration of the latter’s capability in mimicking their biological counterpart.
In addition, for underwater applications, localization of dipole sources has implications in detecting and estimating nearby ﬁsh-like robots and therefore is relevant to robot coordination
and control.
Dagamseh et al. proposed the use of characteristic points (zero-crossings, maxima, etc.)
in the measured velocity proﬁle along the lateral line for dipole source localization [27],
similar to what was proposed by Franosch et al. for modeling the localization by the clawed
4

frog Xenopus [34]. However, this approach would require prohibitively many sensors to
determine the characteristic points, and it is limited to a maximum detection distance of
√
1/ 2 body length (BL). Franosch et al. also suggested a maximum-likelihood estimator-type
model for dipole localization by Xenopus [34], although no algorithms were presented to solve
for the optimal source location. Data-matching/table-lookup approaches were presented by
Pandya et al., where the measured signal pattern was compared with a large, pre-obtained
set of templates or an empirical model ﬁtted with suﬃcient amount of data [68]. These
approaches suﬀered from the need for excessive computing and storage resources, or the
diﬃculty in system-level implementation [68]. Recently, a beamforming algorithm for array
signal processing was used to localize a dipole source and a tail-ﬂicking crayﬁsh, and showed
better performance than the matched ﬁlter [98].

1.2

Research Objective

Inspired by the sensing capabilities of biological lateral lines, our goal in this dissertation
is to investigate signal processing approaches for an artiﬁcial lateral line, which will be
ultimately mounted on an underwater robot such as a robotic ﬁsh (Fig. 1.2). Realization of
this goal could have signiﬁcant implications. First, it will directly beneﬁt many underwater
applications. For example, autonomous underwater vehicles can navigate in narrow spaces
with a near-ﬁeld sensing capability. Second, by learning the mechanisms used by most ﬁsh
and underwater species, we will be able to develop new algorithms for sensor signal processing
that will ﬁnd use beyond hydrodynamic sensing.
5

Figure 1.2: Illustration ot the potential use of an artiﬁcial lateral line system on a robotic
ﬁsh.

1.3
1.3.1

Contributions
An Artiﬁcial Lateral Line System Using IPMC Sensor Arrays

The ﬁrst contribution of this dissertation is a novel approach to the realization of artiﬁcial
lateral lines, which exploits the inherent sensing properties of ionic polymer-metal composites (IPMCs)[81, 19, 58]. As illustrated in Fig. 1.3, An IPMC has three layers, with an
ion-exchange polymer membrane sandwiched by metal electrodes. Inside the polymer, (negatively charged) anions covalently ﬁxed to polymer chains are balanced by mobile, (positively
charged) cations. Deformation under a mechanical perturbation redistributes the cations,
producing a detectable electric signal (e.g., short-circuit current) that is well correlated
with the mechanical stimulus, which explains the sensing principle of an IPMC. Conversely,
under an applied voltage, transport of hydrated cations and water molecules within the
membrane, together with the associated electrostatic interactions, leads to bending motion
6

of the IPMC sample. Many researchers have studied the fabrication [49, 6, 18], characterization, and modeling [63, 86, 10, 17, 14, 28, 71, 72, 73, 100] of IPMC sensors and actuators.
Recent years have seen signiﬁcant interest in using IPMC materials for underwater actuation
[87, 15, 7, 99, 12, 48, 88], sensing [35], and energy harvesting [8].
Under deformation:

Without deformation:
Metal layers
on surfaces

+

+
_

_

_

+

_

+

_

+

_

+

+

_

_

+
+

_

+

_

_

_

+

+

Polymer layer
_

fixed anion

+

+

mobile cation

water

hydrated cation

Figure 1.3: Illustration of IPMC sensing principle.

We have chosen the IPMC material for creating artiﬁcial lateral line in this work for
several reasons. First, IPMC works well in water and has direct mechanosensory property,
which makes it relatively easy to construct the sensor and its readout circuitry. Another
advantage of IPMC sensors (over, e.g., hot-wire ﬂow sensors) is that they automatically
capture the ﬂow polarity. Furthermore, the softness of IPMC material allows it to respond
to small ﬂows and thus attain high measurement sensitivity. Finally, the fabrication processes
for IPMCs are relatively simple and have been studied extensively over the past decade. In
particular, for millimeter-scale and above, one can fabricate an IPMC sheet using well known
7

processes [49] and then cut into any desired sizes.

1.3.2

Localization of a Fixed Dipole Using Neural Network

In this dissertation, we propose an artiﬁcial neural network-based scheme for processing the
signals from IPMC sensor arrays, to localize a dipole source. In general, the relationship
between the source location and the resulting ﬂow ﬁeld (and thus the sensor outputs) is
highly nonlinear and complex. While there is an analytical but nonlinear model for a dipolegenerated ﬂow ﬁeld under ideal assumptions [53], this model may not capture the ﬂow ﬁeld
well under more realistic conditions, an example of which is the frequent scenario that involves ﬂuid-structure interactions. Artiﬁcial neural networks are capable of approximating
nonlinear mappings, and in particular, when such mappings are unknown, there are established methods for training the approximating neural networks [79]. For this reason, we have
adopted a neural network-based scheme for the dipole-source localization. The concept of an
IPMC lateral line and its neural network-based processing scheme are examined experimentally. A prototype comprising six millimeter-scale IPMC sensors, with a body length (BL) of
10 cm, is constructed. The localization experiments are conducted for a dipole source that
is vibrating at 40 Hz. Note that this frequency is consistent with the typical range of dipole
frequencies adopted in the study of biological and artiﬁcial lateral lines (e.g., 50 Hz in [24]
and 45 Hz in [98, 65]).

1.3.3

Nonlinear Estimation-based Dipole Source Localization

An engineering equivalent of a biological lateral line is of great interest to the navigation and
control of underwater robots and vehicles. A vibrating sphere, also known as a dipole source,
8

can emulate the rhythmic movement of ﬁns and body appendages, and has been widely used
as a stimulus in the study of biological lateral lines. Dipole source localization has also become a benchmark problem in the development of artiﬁcial lateral lines. In this dissertation
we present two novel iterative schemes, called Gauss-Newton (GN) and Newton-Raphson
(NR) algorithms, for simultaneously localizing a dipole source and estimating its vibration
amplitude and orientation, based on the analytical model for a dipole-generated ﬂow ﬁeld.
The performance of the GN and NR methods is ﬁrst conﬁrmed with simulation results and
the Cramer-Rao bound (CRB) analysis. Experiments are further conducted on an artiﬁcial lateral line prototype, consisting of six millimeter-scale ionic polymer-metal composite
(IPMC) sensors with intra-sensor spacing optimized with CRB analysis. Consistent with
simulation results, the experimental results show that both GN and NR schemes are able to
simultaneously estimate the source location, vibration amplitude and orientation with comparable precision. Speciﬁcally, the maximum localization error is less than 5% of the body
length (BL) when the source is within the distance of one BL. Experimental results have
also shown that the proposed schemes are superior to the beamforming method, one of the
most competitive approaches reported in literature, in terms of accuracy and computational
eﬃciency.

1.3.4

Localization of a Moving Dipole Source Underwater

We further investigate the problem of localizing a moving dipole source using an artiﬁcial
lateral line (ALL). Such a moving source could represent a swimming ﬁsh or robotic ﬁsh,
which demonstrates a combination of translational motion and oscillatory body/ﬁn motion.
First, we formulate a nonlinear estimation problem based on an analytical model for the
9

moving dipole-generated ﬂow ﬁeld, where we assume that the source location, vibration
amplitude, and moving speed are unknown. We also assume that the amplitudes of the
ﬂow velocities at the DC and dipole vibration frequencies are available as measurements.
Such amplitudes can be obtained by conducting sliding digital Fourier transform (SDFT) of
the measured signal in the time domain. The estimation problem is solved recursively with
the Gauss-Newton method, and we illustrate the proposed approach with simulation and
experimental results.

1.3.5

Underwater Tracking and Size-Estimation of a Moving Object

We investigate the problem of tracking a moving but non-vibrating cylindrical object and
estimating its size and shape using an artiﬁcial lateral line system. Instead of using pressure
sensing as done in the literature for similar problems [11, 33, 91], we consider the measurement of ﬂow velocities as is believed to be what is adopted in biological lateral lines. Based
on a nonlinear analytical model for the moving object-induced ﬂow ﬁeld, two-stage nonlinear
ﬁltering algorithms are proposed to estimate the location, velocity, size, and shape of the
object. Using two-stage ﬁltering instead of one stage is motivated partly by the progressive
information extraction conjectured for biological lateral lines [11], and partly by computational eﬃciency. The approach is ﬁrst illustrated with simulation results, where a moving
cylinder with ellipsoidal cross-section is tracked and its shape parameter estimated. Experimental validation is then conducted with an artiﬁcial lateral line prototype comprising six
IPMC ﬂow sensors. Experimental results show that the proposed approach is able to track
the movement of a cylindrical and an ellipsoidal objects with a circular proﬁle and provide
10

an estimate of the object size, and shape, respectively.

1.3.6

Localization and Identiﬁcation of an Stationary Obstacle
Underwater

Fish exploits the vortex streets generated by leading ﬁsh (schooling behaviors) or by a
static obstacle in ﬂow. First, a ﬁsh swimming between the vortex streets created by two
leading ﬁsh will experience a ﬂow assisting its forward motion (Fig. 1.4(a)) [95]. Second, a
ﬁsh can extracting energy from oncoming vortices. For example, if the caudal ﬁn counterspins a vortex shed by leading ﬁsh or a static obstacle, the energy wasted in the wake can
be minimized (Fig. 1.4(b)) [90]. Such a strategy for autonomous underwater vehicles will
require detecting, tracking and manipulating oncoming vortices using artiﬁcial lateral lines.

For a robot-like ﬁsh or an underwater vehicle, the information sensed via the artiﬁcial
lateral line reﬂects the combined eﬀect of an unknown number of vortices, the ﬂow due to
its own motions, and other disturbances. Biologically, ﬁsh are able to easily identify stimuli
of interest in a noisy environment through its lateral line system [20]. We experimentally
explore the use of IPMC sensors for estimating the hydrodynamic parameters involved in a
K´rm´n vortex street that is created by a stationary cylinder in a ﬂow. We validate that
a a
the vortex shedding frequency, which can be extracted from the sensor signal, shows clear
correlation with the ﬂow speed and the obstacle size.
11

Oncoming vortex
Leading fish

Static obstacle

(a)

Oncoming vortex

Fin-shed vortex
=

+

(b)
Figure 1.4: Beneﬁt of exploiting vortex streets (adapted from [57]): (a) Exploiting ﬂow due
to vortex streets of leading ﬁsh; (b) harvesting energy from oncoming vortices.

1.4

Organization of This Dissertation

The remainder of this dissertation is organized as follows. Chapter 2 presents a complete
development of an IPMC artiﬁcial lateral line. Using IPMC sensors designed and fabricated
12

by our Smart Microsystems Laboratory in Michigan State University (MSU), we assemble
and model the response of an IPMC sensor array, and conduct the sensor calibration. We
then perform fast Fourier transform to analyze the signal spectrum and/or extract signal
amplitudes at frequencies of interest. The performance of the system is demonstrated and
tested through a range of experiments involving various hydrodynamic stimuli. We further
demonstrate the performance of the IPMC lateral line in localizing a dipole source, where an
artiﬁcial neural network is used for signal processing. Chapter 3 ﬁrst reviews an analytical
model for a dipole-generated ﬂow ﬁeld. The digital particle image velocimetry (DPIV) system
(LaVision, Ypsilanti, MI) and Cramer-Rao bound optimization technique are also introduced
in this chapter. New nonlinear estimation algorithms are investigated for localizing a dipole
source with unknown vibration amplitude and orientation. Both simulation and experimental
are presented to support the proposed methods. Chapter 4 and Chapter 5 focus on tracking of
a moving dipole source and a moving but non-vibrating object, respectively. For the moving
dipole case, we demonstrate in both simulation and experimentation that the Gauss-Newton
scheme is able to successfully accomplish the localization and tracking task. For a moving
but non-vibrating cylinder, we propose several two-stage Kalman-ﬁltering-type algorithms,
and conduct simulation and experiments to examine their performance. Chapter 6 presents
the experimental results on how an IPMC lateral line might be used for hydrodynamic
imaging of vortex streets shed by a stationary object in a ﬂow. Chapter 7 draws conclusions
on the research, while the recommendations for future work on the project are proposed in
Chapter 8.

13

Chapter 2
Dipole Source Localization Using
Neural Networks
The signals from IPMC sensors of the artiﬁcial lateral line are complex functions of the
source location. Even for the dipole stimulation, the ﬂow ﬁeld measured with digital particle
image velocimetry (DPIV) deviates appreciably from the theoretical predictions, because
of nonideal ﬂuid conditions and interactions of ﬂuid with structures (e.g., with the IPMC
beams and the walls of the tank). In addition, the sensing characteristics of individual sensors
could be diﬀerent from each other in practice because of imperfect fabrication processes. The
sensor outputs are further contaminated with noises due to ambient water movement and
thermal ﬂuctuations [35]. As a result, it is diﬃcult to decode the sensor signals analytically.
Biological ﬁsh are faced with similar challenges in extracting relevant sensing information
from vast amount of data that are corrupted by noises. However, they manage to accomplish
source localization and other missions robustly through neural network-based information
processing. Taking this biological inspiration, we construct an (artiﬁcial) neural network to
14

process the signals acquired by the IPMC-based lateral line.
The remainder of the chapter is organized as follows. The model for an IPMC sensor
interacting with a ﬂow is ﬁrst discussed in Section 2.1. The neural network-based processing
scheme and the optimal design of the neural network using genetic algorithm are described
in Section 2.2 and Section 2.3, respectively. The experimental setup and basic sensor characterization are presented in Section 2.4. Experimental results on the localization of a dipole
source are shown in Section 2.5. Finally, concluding remarks are provided in Section 2.6.

2.1

Hydrodynamic Model for an IPMC Sensor Under
Oscillatory Flow

In a biological lateral line system, under an AC ﬂow, the ﬂow ﬁeld around the cilium is very
diﬀerent from what it would be under a DC ﬂow. Due to a viscous ﬂow, a boundary layer
symmetric about the cilium will develop at the substrate (Fig. 2.1), and the drag force acting
on the cilium will have some phase shift along the length.

As the starting point of this

analysis, we assume that there is an oscillatory ﬂow with angular frequency ω, and the local
ﬂuid velocity at the ith sensor site is denoted as v(xi , yi ), where sensor i has dimensions of
h × W × t (height, width, thickness) and is located at (xi , yi ). Under the small-deﬂection
assumption, the Euler-Bernoulli equation for the IPMC beam is given as,

YJ

∂ 4q
∂ 2q
+ ρm 2 + c 1
¯
∂y 4
∂t

∂q
−U
∂t
15

+ ρa
¯

∂ 2 q ∂U
−
∂t
∂t2

= ρw
¯

∂U
,
∂t

(2.1)

Figure 2.1: Illustration of the beam dynamics of an IPMC sensor in ﬂow.
where the ﬁrst term represents elastic force (Fe ) acting on the beam [45]:

Fe = Y J

∂ 4q
,
∂y 4

(2.2)

Y is the Yong’s modulus of the IPMC beam, J is the second moment of inertia J = 2 W t3 .
3
The second term in (2.1) represents the inertial force Fm , given as [61]

F m = ρm
¯

∂ 2q
,
∂t2

(2.3)

where ρm = ρm W t, and ρm is the density of the IPMC material. The model treats the
¯
stimulus as an oscillating pressure ﬁeld and we can calculate the ﬂow near the boundary
16

over the ﬂat surface as [38]
−y

U = v(xi , yi ) 1 − e δ

−y
δ

cos

,

(2.4)

where δ is the boundary layer thickness. This is calculated by [38]

2µ
,
ρw ω

δ=

(2.5)

where ρw and µ are the density and dynamic viscosity of water, respectively. The viscous
drag force Fu , on an element of the IPMC beam, is given as the third term in (2.1) [61]

F u = c1

∂q
−U
∂t

,

(2.6)

where c1 = 4πµk, and k is a viscous drag coeﬃcient, which is calculated as

k=

L
L2 + (π/4)2

,

(2.7)

where L is the hydrodynamic force coeﬃcient. L can be calculated as

L = γ + ln

W
√
2δ

,

(2.8)

where γ is Euler’s constant (γ ≈ 0.5772). The acceleration reaction force Fa , is represented
as the fourth term in (2.1) [37]

F a = ρa
¯

∂ 2 q ∂U
−
∂t
∂t2

17

.

(2.9)

Pressure ﬁeld that develops around the IPMC beam as ﬂow accelerates results in a buoyant
force, Fb , that can be calculated as

F b = ρw
¯

∂U
,
∂t

(2.10)

where ρw = ρw W t. By converting the equation (2.1) from the time-domain to the Laplace
¯
(s) domain, we can derive the transfer function that relates the beam tip deformation q1
to the local velocity v. With typical parameters, the transfer function turns out to be
approximated well by
q1 (s)
ca
= ,
v(s)
s

(2.11)

where ca is constant. On the other hand, Chen et al. [16] have developed a dynamic for
IPMC sensors, where the transfer function from the tip displacement q1 to IPMC current
output I can be approximated as
I(s)
= cb s,
q1 (s)

(2.12)

where cb is a constant. The overall approximate model from local velocity v to the IPMC
current output I is then given as
I
= η,
v(s)

(2.13)

where η = ca cb is constant gain that will be calibrated in experiments. Note that the approximations used to arrive at (2.11) and (2.12) are only valid for an intermediate frequency
range, where the IPMC sensor behaves like a band-pass ﬁlter. For DC or high-frequency
range, the sensor model will need to be adapted appropriately.

18

2.2

Source Localization Using Neural Network Processing

We consider the localization of a dipole source in a two-dimensional (2D) plane with coordinates (x, y). The proposed approach can be extended readily to the localization in a
3D space, although that will require more training points for the neural network. Fig. 2.2
illustrates the problem setup. Here (x, y) denotes the location of the source, and the IPMC
sensor arrays are placed at a known location.

As illustrated in Fig. 2.3, we adopt the

y
Dipole source
at ( x, y )

IPMC sensor array

x
Figure 2.2: Schematic of the dipole source localization problem.

multilayer perceptron (MLP) architecture for the neural network. An MLP network consists
of an input layer, a hidden layer, and an output layer, and is the most widely used network
structure for nonlinear classiﬁcation and prediction applications [79].
One could use diﬀerent features extracted from the sensor output data as the input to
the neural network. In this work, we use the signal amplitude at the stimulus frequency
19

Sensor 1

x

Sensor 2

y
Sensor N

Output layer
Input layer
Hidden layer

Figure 2.3: Schematic of the MLP neural network for signal processing of the IPMC lateral
line.
because of its robustness to measurement noises. The amplitude is obtained through fast
Fourier transform (FFT) of the raw sensor signal. The number of inputs is the same as the
number of IPMC sensors considered. For comparison purposes, in this work we investigate
the performance of the artiﬁcial lateral line when diﬀerent numbers of sensors are adopted.
The number of the hidden-layer nodes is chosen through a genetic algorithm (GA)-based
optimization process, which will be further described below. Each hidden layer node represents the operation of nonlinear activation, which takes the form of a sigmoid function. The
output layer has two nodes, representing the x and y coordinates of the vibrating source.
The number of hidden-layer nodes and the connective weights between the layers are
determined through a two-phase training procedure, using the software NeuroSolutions [64].
The training data are obtained by placing the stimulus at known locations (xi , yi ), 1 ≤ i ≤
M , measuring the corresponding sensor outputs and computing the signal amplitudes. In
the ﬁrst training phase, a genetic algorithm is used to ﬁnd an appropriate value for the
number of hidden-layer nodes and a reasonable set of values for the connective weights. In
particular, each genome encodes both the number of hidden-layer nodes and the weights of
20

all connecting edges. The maximum number of hidden-layer nodes is limited to 24 based on
the numbers of network inputs and outputs. We consider a population size of 25 genomes
and run 40 generations of evolution (crossover/mutation). The ﬁtness function (or more
appropriately, the cost function) used in the selection process is the mean square error J:

1
J=
2M

M
i=1

(xi − xi )2 + (yi − yi )2 ,
ˆ
ˆ

(2.14)

where (ˆi , yi ) denotes the predicted value for (xi , yi ) under the current network structure
x ˆ
and weights.

The values of the connective weights obtained in the ﬁrst training phase then serve as the
initial condition for weights reﬁnement in the second phase, where the network structure is
ﬁxed as determined in the ﬁrst phase. Delta-bar-delta learning [79], with adaptive learning
rate, is used for weights optimization. Let K be the total number of weights. For each
weight wk , 1 ≤ k ≤ K, the update rule is
new
old
new ∂J ,
w k = w k − ηk
∂wold

(2.15)

k

where the adaptive learning rate ηk is updated as

new
ηk =


 old
 η + a, if
 k






∂J > 0
∂wold
k
,
∂J ≤ 0
old , if
bηk
∂wold
k

and a, b are constants satisfying 0 < a, b < 1.
21

(2.16)

2.3

Optimum Design of the Neural Network Using Genetic Algorithm

The genetic algorithm (GA) is a search technique, based on an abstraction of biological evolution, used to ﬁnd solutions of optimization problems. In GA, a population of chromosomes,
which are represented as binary vectors, evolves into a new population of chromosomes using
the natural selection forces of crossover and mutation. Some chromosomes are selected and
allowed to reproduce and the more “ﬁt” chromosomes produce more oﬀspring. Recombination is achieved via the crossover operator, which exchanges parts of two chromosomes.
Mutation can randomly change some segments in the chromosome. Chromosomes are used
to encode possible solutions to an optimization problem. The quality of a solution is represented by its ﬁtness, which is the objective function to be optimized. In our case, the
objective function is the mean-square error. The GA is implemented according to the following steps: (i) Randomly generate an initial chromosome population; (ii) calculate the
minimum square error to indicate the ﬁtness function value for each chromosome in the
population; (iii) create a new chromosome population from mating of current chromosomes
and applying mutation and recombination; (iv) select a ﬁxed number of chromosomes such
that more ﬁt chromosomes are more likely to be selected; (v) compute the ﬁtness of the new
population; (vi) increment the current number of generation. Finally, if not (end-test) go
to step (iii), or else stop and return the best chromosome. We use the genetic algorithm
to determine the optimal initial weights and the optimum number of nodes in the hidden
layer of the neural network. The optimal initial weights will guarantee a faster convergence
to the global minimum of the objective function, and the optimum number of nodes in the
22

hidden layer will avoid the bias-variance dilemma of neural networks. Given a data set, it is
known that a small number of nodes in the hidden layer captures the data trend only in a
small region of the pattern, whereas a large number of nodes in the hidden layer results in
overﬁtting the data [79].

2.4

Experimental Setup and Sensor Characterization

Fig. 2.4 shows the constructed lateral line prototype, consisting of six IPMC sensors. Each
sensor, with dimensions 8 mm × 2 mm × 200 µm, has been cut from an IPMC sheet
fabricated by the Smart Microsystems Laboratory at Michigan State University, following a
recipe similar to the one described in [49]. The IPMC material contains lithium ions (Li+ )
as cations. The sensor-to-sensor separation is 2 cm, resulting in a total span of 10 cm, which
will be regarded as the Body Length (BL) in later discussions.

Figure 2.4: A prototype of IPMC-based lateral line system.

Under a mechanical stimulus, an open-circuit voltage or a short-circuit current can be
measured across the two electrodes of an IPMC. We have chosen to take the short-circuit current as the sensor output because current measurement is less susceptible to noises. Fig. 2.5
shows the schematic of the measurement circuit, which consists of two cascaded operational
23

ampliﬁers (op-amps)[35]. A low-noise, low-bias precision op-amp (OPA 124 from Texas
Instruments) is adopted for the ﬁrst-tier ampliﬁcation, to reduce both the noise and the
spurious DC bias in the sensor output induced by the leakage current. Since the “−” terminal of Op-amp 1 is virtually the ground (following standard op-amp circuit analysis), the
two electrodes of IPMC can be considered short-circuited. The sensing current generated
under this conﬁguration, i(t), is proportional to the voltage output v1 (t) = R1 i(t). The
second op-amp is introduced for gain adjustment, where the resistor R3 is tunable. The
output v2 (t) is related to the current signal i(t) via v2 (t) =

R3 R1
R2 i(t). In the circuit we

used, R1 = 470 kΩ, R2 = 10 kΩ, and R3 was adjustable from 0 to 50 kΩ. Acquisition and
processing of the IPMC sensor output are conducted through a dSPACE system (DS1104,
dSPACE Inc., Germany). A digital low-pass ﬁlter with cutoﬀ frequency 55 Hz is further
implemented to remove high-frequency noises from the sensor signals .

d (t )
R1

R3

IPMC

i (t )

v1 (t ) R2

v2 (t )

+
+

GND

Op-amp 1
Op-amp 2

Figure 2.5: The schematic of circuit for measuring the short-circuit current generated by the
IPMC sensor.

All experiments are conducted in a water tank, shown in Fig. 2.6(a). The tank measures 6 × 2 × 2 ft3 , and is equipped with a digital particle image velocimetry (DPIV) system
24

(LaVision, Ypsilanti, MI). In a DPIV system, small particles are dispersed in a ﬂuid and a
laser sheet is created in the ﬂuid to illuminate the particles. Processing of images taken in
quick successions can reveal the movement of particles and thus provide information about
the ﬂow ﬁeld. In our experiments, the DPIV system is used for preliminary characterization
of the ﬂow ﬁeld under the stimulating source. A dipole source is created with a mini-shaker
(4810, Br¨el & Kjaer, Denmark) (Fig. 2.6(b)), the vibration amplitude and frequency of
u
which can be readily controlled through a voltage input to the shaker. A lightweight bar
ﬁrmly attached to the mini-shaker then translates the vibration to a sphere rigidly coupled
to the bar. The sphere, which is a steel ball, has a diameter of 19 mm. The frequency range
of the dipole spans DC–20 kHz. The source location and vibration direction with respect to
the IPMC lateral line can be adjusted by moving the stand holding the IPMC lateral line or
by moving the source itself.

Camera

Tank

Laser
(a)

(b)

Figure 2.6: Experimental setup: (a) Experimental tank equipped with a DPIV system; (b)
the mini-shaker-based dipole source.

DPIV experiments have been conducted to observe the ﬂow ﬁeld around the IPMC
sensors, for diﬀerent source-sensor separation, to get a sense of the ﬂow magnitude un25

der the stimuli. For example, Fig. 2.7 shows the ﬂow ﬁelds around an IPMC sensor for a
19 mm × 14 mm window, under a 20 Hz stimulus produced with a dipole source described
in [35], where both the source vibration axis and the illuminated plane are perpendicular to
the IPMC beam plane. In other words, the IPMC beam would bend left and right in Fig. 2.7,
although the IPMC thickness did not look to scale because of laser light scattering around
the beam. Fig. 2.7(a) and (b) correspond to dipole-IPMC separations of 10 and 45 cm,
respectively, which clearly indicates that that the ﬂow velocity decreases as the separation
between source and sensor increases.

(a)

(b)

Figure 2.7: Flow ﬁeld around one IPMC sensor under a dipole stimulation: (a) Dipole-sensor
separation 10 cm, with maximum ﬂow amplitude of 9.3 mm/s in the ﬁeld; (b) dipole-sensor
separation 45 cm, with maximum ﬂow amplitude of 4.1 mm/s in the ﬁeld.

Fig. 2.8 shows typical sensor responses under the mini-shaker-based stimulation (40Hz),
indicating that the current output from the IPMC sensor is at the order of µA, which can
be captured very well by the sensing circuit. Fig. 2.8(a) shows the signal from sensor #2
of the IPMC lateral line, when the source is 1 cm away from this sensor, while Fig. 2.8(b)
shows the response from sensor #4 (4 cm from sensor #2) under the same stimulus. In this
26

experiment the lateral line axis is perpendicular to the vibration axis of the dipole source,

Sensor output (µ A)

and perpendicular to the beam plane of each sensor.

2
1
0
−1
−2
9

9.1

9.2

9.3

Time (s)

9.4

9.5

Sensor output (µ A)

(a)
0.3
0.2
0.1
0

−0.1
−0.2
9

9.1

9.2

9.3

9.4

9.5

Time (s)

(b)
Figure 2.8: Typical IPMC sensor signals: (a) Signal from sensor #2, which is close to the
dipole source; (b) signal from sensor #4.

Fig. 2.9 shows the amplitude of the signal from one IPMC, located at (7,0), as a function
of the location (x, y) of the mini-shaker-based dipole stimulus vibrating parallel and perpendicular to the IPMC beam plane.

It can be seen that, while in general the signal gets

stronger when the source gets closer, the overall amplitude landscape has a sophisticated
proﬁle. The results in Fig. 2.9 clearly illustrate the challenge in underwater localization.
In particular, it would be diﬃcult to localize unambiguously a source with a single sensor;
27

Sensor output (µ A)

1
0.8
0.6
0.4
0.2
0
10

Sensor
20
15

5

10

y (cm)

x (cm)

5
0

0

(a)

Sensor output (µ A)

1
0.8
0.6
0.4
0.2
0
10

Sensor
20
15

5

y (cm)

10
0 0

5

x (cm)

(b)
Figure 2.9: Sensor characterization, measured sensor signal magnitude as a function of source
location, The dipole vibration axis is: (a) Parallel to the IPMC beam plane; (b) perpendicular
to the IPMC beam plane.

28

instead, a lateral line-like array structure will be needed.

2.5

Experimental Results on Diploe Source Localization

As shown in Fig. 2.11, a working area of about 20 × 10 cm2 is used in the experiments.
The frequency of the dipole source is 40 Hz, the vibration is along the y direction with
an amplitude of 1.91 mm. The lateral line is placed in parallel with the x-axis, at y = 0
and spanning from x = 5 cm to x = 15 cm. The beam plane of each IPMC sensor is
perpendicular to the lateral line direction. Both the dipole source and the lateral line are
placed 6 cm underwater. First, for the purpose of training the neural network, we deﬁne a grid
of 1 × 1 cm2 lattices in the working area, resulting in a total of 180 grid points. The dipole
source is placed at each grid point, where the amplitude of the signal from each IPMC sensor
in the lateral line is obtained. Subsets of these data are used to train the neural network in
diﬀerent cases. Fig. 2.10 shows the evolution of the cost function J in (2.14) during the ﬁrst,
genetic algorithm-based training phase, for the case involving six sensors and 180 training
points. It can be seen that the cost function converges to a steady-state value as the number
of generations approaches 40. As explained in Section 2.2, the connective weights obtained
in the ﬁrst training phase are then reﬁned in the second phase using delta-bar-delta learning.
To test the performance of the lateral line, we place the dipole source at each of 16 points
along an elliptical track centered at (10, 5.5), obtain the amplitudes of sensor signals from
the lateral line and feed them to the trained neural network, and compare the output of the
network to the actual location.
29

1

Fitness value

0.8
0.6
0.4
0.2
0
0

10

20

30

40

Number of generations

Figure 2.10: Evolution of the ﬁtness function (the mean square error J) over generations,
in the genetic algorithm-based training phase. This example involves all six sensors and 180
training points.
Fig. 2.11(a) and (b) show the localization results when all six sensors are used, where
90 (2 × 1 cm2 lattices) and 50 (2 × 2 cm2 lattices) training points are utilized, respectively.
The locations of the training points are denoted with squares in the ﬁgures. Fig. 2.12(a) and
(b) further show the actual localization errors at all 16 validation points for the two cases.
It is clear from Fig. 2.12(b) that, for a training grid of 2 × 2 cm2 lattices (50 points), the
maximum localization error is 0.3 cm. The maximum error drops to 0.2 cm when the number
of training points is increased to 90 (Fig. 2.12(a)). Note that the resulting neural networks
are diﬀerent when diﬀerent numbers of training points (or diﬀerent numbers of sensors) are
used. For example, the neural network has 13 hidden nodes for the case in Fig. 2.11(a),
while it has 10 hidden nodes for the case in Fig. 2.11(b).
In order to examine the eﬀect of the number of sensors on the localization performance,
we have repeated the localization with four and two sensors of the lateral line, as shown
30

Validation: actual
Validation: predicated
Training
10

y (cm)

8
6
4
2
0
0

5

Sensors
10

15

20

x (cm)

(a)
Validation: actual
Validation: predicated
Training

10

y (cm)

8
6
4
2
0
0

5

Sensors
10

15

20

x (cm)

(b)
Figure 2.11: Localization of the dipole source using all six sensors of the lateral line: (a)
With 90 training points; (b) with 50 training points.

in Figs. 2.13 and 2.14, respectively. In each case, the sensors are taken from the middle
section of the lateral line. It is clear that with fewer sensors, the localization error increases.
Fig. 2.15 further shows the maximum and average localization errors when diﬀerent numbers
31

Error (cm)

0.2
0.15
0.1
0.05
0
0

2

4

6

8

10

12

14

16

Index of points

(a)

Error (cm)

0.3
0.2
0.1
0
0

2

4

6

8

10

12

14

16

Index of points

(b)
Figure 2.12: Localization error at 16 validation points when all six sensors are used: (a)
With 90 training points; (b) with 50 training points.

of sensors and training points are used. From the ﬁgure, we can see that the localization
performance with four sensors becomes comparable to that with six sensors when the number
of training points is quadrupled (from 2 × 2 cm2 lattices to 1 × 1 cm2 lattices). And for two
sensors, even with all 180 training points, the maximum localization error is bigger than that
with all six sensors when only 25 training points are used.

2.6

Chapter Summary

In this chapter, we investigated a new approach to the realization of artiﬁcial lateral lines
for underwater robots and vehicles, including both the proposal of using IPMC materials as
32

sensing elements and the neural network-based signal processing algorithm. The eﬀectiveness
of the proposed approach was validated in experiments involving localization of a dipole
source. Experimental results show that the IPMC-based lateral line can localize the dipole
source 1-2 BLs away, with a maximum localization error of 0.3 cm, when the data for training
the neural network are collected from a grid of 2 cm by 2 cm lattices. The performance of the
lateral line using fewer sensors is also studied, and it is found that, with fewer sensors, more
training data are needed to achieve a given localization precision. For example, when four
sensors are used, it requires training data on a grid of 1 cm by 1 cm lattices to yield the same
level of localization accuracy as achieved by six sensors with training data from 2 × 2 cm2
lattices. These results have not only demonstrated the feasibility of using arrays of IPMC
sensors as artiﬁcial lateral lines, but also provided interesting insight into the tradeoﬀs in
design and implementation.

33

Validation: actual
Validation: predicated
Training
10

y (cm)

8
6
4
2

Sensors

0
0

5

10

15

20

x (cm)
(a)
Validation: actual
Validation: predicated
Training
10

y (cm)

8
6
4
2
0
0

Sensors
5

10

15

20

x (cm)
(b)
Figure 2.13: Localization of the dipole source using four sensors of the lateral line: (a) With
90 training points; (b) with 50 training points.
34

Validation: actual
Validation: predicated
Training
10

y (cm)

8
6
4
2
Sensors

0
0

5

10

x (cm)

15

20

(a)
Validation: actual
Validation: predicated
Training
10

y (cm)

8
6
4
2
0
0

Sensors
5

10 x (cm)

15

20

(b)
Figure 2.14: Localization of the dipole source using two sensors of the lateral line: (a) With
90 training points; (b) with 50 training points.

35

Max prediction error (cm)

6 Sensors
4 Sensors
2 Sensors

1.5

1

0.5

0
0

50

100

150

200

Number of training points

Average prediction error (cm)

(a)
6 Sensors
4 Sensors
2 Sensors

0.6
0.5
0.4
0.3
0.2
0.1
0
0

50

100

150

200

Number of training points
(b)
Figure 2.15: Error in localization of the dipole source using the lateral line: (a) Maximum
error along the track; (b) average error along the track.

36

Chapter 3
Nonlinear Estimation-based Dipole
Source Localization
Most existing work deals with a dipole source that has ﬁxed vibration amplitude and orientation. In this chapter we consider the dipole localization problem where both the vibration
amplitude and vibration orientation are unknown. When the vibration amplitude is varying,
the problem becomes challenging, partly because a source far away but with large vibration
and another source nearby but with small vibration could produce signals of similar amplitudes. If we consider the scenario of one ﬁsh detecting a predator using its lateral line system,
the solution to the aforementioned problem would imply that the ﬁsh not only knows where
the predator is, but also which way it is swimming toward. There are similar implications
for artiﬁcial lateral line systems.
To solve these problems, we treat the source location, vibration amplitude, and orientation of vibration axis as unknowns of a nonlinear estimation problem, and then present two
recursive algorithms. In the ﬁrst scheme, we solve a linearized estimation problem with the
37

Gauss-Newton (GN) method. In the second scheme, we solve the nonlinear equations corresponding to the ﬁrst-order optimality condition using the Newton-Raphson (NR) method.
We will also compare the proposed methods with the beamforming method proposed by [98],
which is among the most competitive methods reported in the literature. In this chapter we
also exploit the Cramer-Rao Bound (CRB) for the assessment of estimation performance.
The remainder of the chapter is organized as follows. The estimation problem is formulated in Section 3.1. The GN, NR, and beamforming algorithms are described in Section 3.2, along with the CRB analysis. Simulation and experimental results are presented
in Section 3.3 and Section 3.4, respectively. Finally, chapter conclusions are provided in
Section 3.5.

3.1

Problem Formulation

As widely assumed in the literature [46, 39, 43], we consider a potential ﬂow generated by
a vibrating sphere, the dipole source. The velocity potential ϕ at a spatial point can be
expressed as [53]
ϕ(r) =

a3 (vd · r)
,
2 r 3

(3.1)

where a is the diameter of vibrating sphere, vd is the instantaneous velocity of the dipole
source, r represents the relative location of the point of interest with respect to the dipole,
and

·

denotes the Euclidean norm of a vector. At the point r, the ﬂow velocity v(r) is

v(r) = −∇ϕ(r) =

a3 3(vd · r)r − r 2 vd
.
2 r 5
38

(3.2)

Since our goal in this dissertation is to localize the dipole source, we redeﬁne r hereafter as
the location of the dipole source with respect to a known point. Note that the formula (3.2)
for the dipole-generated velocity ﬁeld remains valid under the new deﬁnition of r, because
v(−r) = v(r).
To ease the discussion, we will focus on the case where the lateral line sensor and the
dipole source are located in the same plane (denoted as the x − y plane) in which the axis of
dipole oscillation lies; see Fig. 3.1 for illustration. This is a scenario widely adopted in the
study of biological and artiﬁcial lateral lines, and our proposed method will extend to the
more general, three-dimensional case in a straightforward manner. As depicted in Fig. 3.1,
the dipole source is located at (xs , ys ) and vibrates with velocity v d , where the angle between
the vibration axis and the x−axis is φ. Let an array of N sensors be located in parallel to
the x-axis, with their locations assumed to be known and denoted as (xi , yi ), 0 ≤ i ≤ N − 1.
We assume that the presence of sensors has negligible eﬀect on the ﬂow distribution as
characterized in (3.2); further discussion on sensor-ﬂow interactions will be provided in
Subsection 3.4.3. We assume that each sensor is able to provide a noisy measurement of
the local ﬂow velocity v along one direction, which, without loss of generality, is taken to be
the x−direction in this dissertation. The latter consideration is motivated by the fact that
most beam-like ﬂow sensors [92, 74, 31, 60, 97, 78, 29, 1] can only sense one component of
the ﬂow velocity. In the absence of noise, from (3.2), we get

v(xi , yi ) =

a3 v d
((2(xi − xs )2 − (yi − ys )2 ) cos φ
5
2 ri
+ 3(xi − xs )(yi − ys ) sin φ),

where ri = (xi − xs , yi − ys )T , with “T ” denoting the transpose.
39

(3.3)

Let the sphere vibrate

y
( xs , ys )

Dipole source

ri

( xi , yi )
Sensor i

x
Figure 3.1: Illustration of the problem setup. The dipole source is located at (xs , ys ).

with angular frequency ω and amplitude A, i.e.,

v d = A sin(ωt).

Furthermore, due to the periodic nature of the dipole-generated ﬂow, we consider extracting
the signal amplitude at frequency ω from each sensor as the measurement. In particular, given
n

2
the raw signals from time n1 to n2 , {si (n)}n=n1 , the signal amplitude at ω can be readily

obtained through techniques such as fast Fourier transform (FFT) [66] or sliding discrete
Fourier transform (SDFT) [44, 82], and this inherent ﬁltering process greatly reduces the
noise contained in the raw sensor signal. For sensor i, the measurement Mi can then be
40

written as

Mi =

a3 ((2(xi − xs )2 − (yi − ys )2 )α1 + 3(xi − xs )(yi − ys )α2 )
+ di ,
2 ri 5

(3.4)

where α1 = A cos φ, α2 = A sin φ, and di is the (amplitude) measurement noise for sensor
i. Note that, to accommodate the nonnegativity constraint for the extracted amplitude, we
have placed the noise di inside | · |, since in general the noise could be positive or negative.
N −1
The noises {di }i=0 are assumed to have zero means and be uncorrelated with each other.

Based on the sensor measurements, we are interested in estimating the location (xs , ys ),
vibration amplitude A, and vibration direction φ of the dipole source. For ease of presentation, for i = 0, · · · , N − 1, deﬁne
fi (θ) =

a3
((2(xi − xs )2 − (yi − ys )2 )α1 + 3(xi − xs )(yi − ys )α2 ),
2 ri 5

(3.5)

where θ = (α1 , α2 , xs , ys ) represents the set of parameters to be estimated. Once θ is
obtained, we can ﬁnd the vibration amplitude A =

2
2
α1 + α2 and vibration orientation

φ = arctan(α2 /α1 )). We further deﬁne

△

M = (M0 , · · · , MN −1 )T ,
△

f (θ) = (f0 (θ), · · · , fN −1 (θ))T ,
△

d = (d0 , · · · , dN −1 )T .
We then have
M = |f (θ) + d| ,
41

(3.6)

where | · | denotes the component-wise absolute value of the vector. The estimation problem
is formulated as follows. Let the dipole diameter a, frequency ω, and the sensor locations
N −1
ˆ
{(xi , yi )}i=0 be known. Given the sensor measurement M , provide an estimate θ for the

parameter θ, such that

ˆ
ˆ
J(θ) = M − f (θ)

T

ˆ
M − f (θ)

(3.7)

ˆ
is minimized, where f (θ) is the predicted measurement based on the estimated parameters.
Before presenting the estimation algorithms, we brieﬂy comment on the assumptions
N −1
made on the knowledge of a, ω, and {(xi , yi )}i=0 . Clearly, it is natural to assume knowing
N −1
{(xi , yi )}i=0 . If the frequency ω is not given a priori, it can be identiﬁed through FFT

of the sensor signals and scan to capture the frequency that has a maximum magnitude.
Finally, if a is unknown, one will be able to estimate a3 A as a lump parameter, since that is
how a and A collectively aﬀect the ﬂow ﬁeld.

3.2

Estimation Algorithms

In this dissertation we propose and compare two recursive algorithms to obtain the parameˆ
ter estimate θ for minimizing the cost function J in (3.7). The function f (θ) is linear in α1
ˆ
and α2 , but highly nonlinear in (ˆs , ys ). In the ﬁrst algorithm, brieﬂy summarized in Subsecˆ
x ˆ
tion 3.2.1, |f (θ)| is linearized at each intermediate parameter estimate, and the estimate at
the next iteration is obtained by solving a standard least-squares problem. This algorithm
is also known as the Gauss-Newton method [47]. In the second algorithm, described in
Subsection 3.2.2, we ﬁrst derive the (nonlinear) equation capturing the ﬁrst-order necessary
42

condition for minimizing J, and then solve the equation using the Newton-Raphson method
[47].
We will further compare our proposed algorithms experimentally with the beamforming
algorithm, which was one of the most competitive algorithms reported for dipole localization
[98]. The algorithm will be brieﬂy reviewed in Subsection 3.2.3.
The Cramer-Rao bound (CRB) [25] establishes the lower limit for the variance of any
unbiased estimator, and we have conducted CRB analysis for both the optimal design of
a lateral line sensor and the evaluation of the proposed algorithms. The CRB analysis is
presented in Subsection 3.2.4.

3.2.1

Gauss-Newton Algorithm

ˆ
¯
Linearizing the function f (θ) about some nominal point θ, we have

△

ˆ
ˆ
¯
¯ ˆ ¯
f (θ) ≈ fL (θ) = f (θ) + B(θ)(θ − θ),
where
¯
B(θ) =

∂ |f (θ)|
.
∂θ θ=θ
¯

ˆ
The estimation problem is then converted to ﬁnding θ, such that

ˆ
ˆ
J1 (θ) = M − fL (θ)
43

T

ˆ
M − fL (θ)

(3.8)

is minimized. This becomes a standard least squares estimation problem, the solution of
which is

ˆ
¯
¯
¯
θ = θ + λ B(θ)T B(θ)

−1

¯
¯
B T (θ) M − f (θ) ),

(3.9)

where λ is a stepping parameter. To minimize the impact of error introduced by the lin¯
earization, a recursive version of (3.9) is used, where θ is replaced by the estimate at the
previous iteration:

ˆ
ˆ
ˆ
ˆ
θk+1 = θk + λ B(θk )T B(θk )

−1

ˆ
ˆ
B T (θk )(M − f (θk ) ),

(3.10)

ˆ
ˆ
with θ1 = θ0 , an initial guess for the parameters. The algorithm is stopped when θk+1 −
ˆ
θk ≤ ǫ, where ǫ > 0 is a speciﬁed tolerance.

3.2.2

Newton-Raphson Algorithm

The necessary condition for minimizing (3.7) is

∂J
= 0,
ˆ
∂θ

(3.11)

which implies
△ ∂ |f | T
ˆ
) (M − f (θ) ) = 0.

ˆ
g(θ) = (

ˆ
∂θ

44

(3.12)

An estimate of the parameter θ can then be obtained by solving the nonlinear equation
(3.12) using the Newton-Raphson method:

ˆ
ˆ
ˆ
ˆ
θk+1 = θk − λG−1 (θk )g(θk ),

(3.13)

ˆ
where λ is stepping parameter, and G(θ) denotes ∂g and can be evaluated as
ˆ
∂θ

ˆ
G(θ) =

N −1
i=0

ˆ
Mi − fi (θ)

∂ 2 |fi |
−
ˆ
∂ θ2

∂ |fi | T
ˆ
∂θ

∂ |fi |
ˆ
∂θ

.

ˆ
The starting value θ1 in (3.13) is set to be θ0 , an initial guess of the parameters. Similar to
ˆ
ˆ
the Gauss-Newton method, the iteration in (3.13) is stopped when θk+1 − θk ≤ ǫ for a
prespeciﬁed ǫ.

3.2.3

Beamforming Algorithm

The beamforming algorithm based on Capon’s method was used by Yang et al. to estimate
a ﬁve-dimensional (5D) dipole state (3D location and 2D vibration orientation) [98]. In
our dipole localization setup, the dipole state is fully captured by the 4D parameter vector
θ. Given the raw sensor signals for instantaneous ﬂow velocities (projected into the sensing
n

2
ˆ
axes), denoted as {s0 (n), · · · , sN −1 (n)}n=n1 , one constructs an energy-like map E(θ), which

ˆ
indicates how likely the actual dipole state is θ. The optimal estimate of θ is then obtained
as
ˆ
ˆ
θ∗ = arg max E(θ).
ˆ
θ

(3.14)

ˆ
In order to evaluate E(θ), one ﬁrst evaluates the correlation matrix R based on the
45

measurements from all sensors:
n

2
1
S(n)S T (n),
R=
n2 − n1 n=n
1

where S(n) = [s0 (n), s1 (n), · · · , sN −1 (n)]T . The energy-like function E is then computed
as
ˆ
E(θ) =

1
ˆ
ˆ
|f (θ)|T R−1 |f (θ)|

.

(3.15)

ˆ
To search for the maximizing θ for E, one typically scans through the space for the dipole
state [98]. The latter involves creating a discrete grid of dipole state components, and the
resulting total number of grid points increases exponentially with the dimension of the dipole
state.

3.2.4

Cramer-Rao Bound Analysis

Consider (3.4). For the Cramer-Rao bound (CRB) analysis, we assume that the noise di
2
from each sensor is Gaussian with zero mean and variance σm . Given the dipole state θ, the

probability for the sensors to have the measurement M = (M0 , · · · , MN −1 )T is expressed as
P (M, θ) =

=

N −1

1

2
(2πσm )1/2
i=0

−1 (M −f (θ))2
−1 (−M −f (θ))2
i i
i i
2
2
2σm
2σm
e
+e




2
2
−Mi −fi (θ)
M fi (θ) 
2
2
2
 2σm 2σm
e
cosh i 2  .
e
2 )1/2
σm
(2πσm
i=0

N −1

46

,

(3.16)

(3.17)

The log-likelihood function is then
N −1

N −1

1
1
Mi2 − 2
fi2 (θ)
ln P (M, θ) = N
− 2
2
2σm
(2πσm )1/2 2σm i=0
i=0
2

+

N −1

ln cosh

i=0

Mi fi (θ)
2
σm

.

(3.18)

Taking the ﬁrst derivative,

∂ ln P (M, θ)
1
= 2
∂θ
σm

N −1
i=0

−fi (θ)

∂fi (θ)
+ tanh
∂θ

Mi fi (θ)
2
σm

Mi

∂fi (θ)
∂θ

,

and consequently,
∂ 2 ln P (M, θ)
∂θ2

=

1
2
σm



N −1
i=0



− fi (θ)

∂ 2 fi (θ)
∂fi (θ)
−
2
∂θ
∂θ

∂fi (θ) T
∂θ

∂ 2 fi (θ) Mi2 ∂fi (θ)
Mi fi (θ)
Mi
+ 2
+ tanh
2
∂θ
σm
∂θ2
σm
M fi (θ)
× 1 − tanh2 i 2
.
σm

47

∂fi (θ) T
∂θ
(3.19)

Taking the negative expectation on (3.19) results in

− E



N −1

∂ 2 ln P (M, θ)
∂θ2

Mi fi (θ)
2
σm
M fi (θ)
1 − tanh2 i 2
σm
N −1
i=0



i=0
∂ 2 fi (θ) Mi2
Mi
+ 2
∂θ2
σm

+ tanh

1
=
2
σm

1
=− 2 E
σm 

1 − tanh2

Mi fi (θ)
2
σm

Mi

Mi fi (θ)
2
σm

∂fi (θ)
∂θ

∂fi (θ)
∂ 2 fi (θ)
−
∂θ
∂θ2

∂fi (θ) T
∂θ

∂fi (θ) T
∂θ

∂fi (θ) T
∂θ

∂ 2 fi (θ)
∂fi (θ)
fi (θ)
+
2
∂θ
∂θ

−E tanh

− fi (θ)

∂fi (θ) T
∂θ

∂ 2 fi (θ) Mi2 ∂fi (θ)
+ 2
∂θ
∂θ2
σm

,

The Fisher information matrix I(θ) can then be derived as
∂ 2 ln P (M ; θ)
∂θ2

N −1
1
∂ 2 fi (θ)
∂fi (θ)
=
fi (θ)
+
2
2

∂θ
σm
∂θ

I(θ) = −E

i=0

−

∞

Mi =0

gi (Mi , θ)Pi (Mi , θ)dMi

∂fi (θ) T
∂θ

,

(3.20)

where E[·] denotes the expectation,

Mi fi (θ)
2
σm

△

gi (Mi , θ) = tanh

M 2 ∂fi (θ)
+ 2i
∂θ
σm
Pi (Mi , θ) =

1
2
(2πσm )1/2

Mi

∂ 2 fi (θ)
∂θ2

M fi (θ)
∂fi (θ) T
1 − tanh2 i 2
,
∂θ
σm

−1 (−M −f (θ))2
−1 (M −f (θ))2
i i
i i
2
2
2σm
2σm
+e
e

48

(3.21)
.

(3.22)

Note that I −1 (θ) represents the lower bound (CRB) on the estimation variance for any
ˆ
ˆ
unbiased estimator of θ; speciﬁcally, the variance of the j-th component of θ, denoted as θ[j] ,
satisﬁes
ˆ
Var(θ[j] ) ≥ I −1 (θ)

j,j

, j = 1, · · · , 4.

(3.23)

In Section 3.3, we will illustrate the use of (3.23), with j = 1, 2, for the optimal design of a
lateral line sensor (with given vibration amplitude and orientation) and for the localization
performance evaluation of the proposed estimation algorithms. Appendix A provides the
ﬁrst and second derivatives of the fi (·).

3.3
3.3.1

Simulation Results
CRB-based Optimal Design of Lateral Line Sensor

The CRB analysis presented in Subsection 3.2.4 can be used for the optimal design of a lateral
line system, including both the number of ﬂow sensors and the geometric arrangement of
these sensors. The latter point is illustrated here with an example of optimizing the intrasensor spacing, which is later adopted in the construction of an experimental prototype. We
consider a working area of 20 × 10 cm2 (−10 ≤ x ≤ 10, 0 < y ≤ 10), and assume that the
lateral line consists of 6 ﬂow sensors, uniformly distributed along the x−axis and centered
at the origin. We would like to determine the sensor-to-sensor spacing that minimizes the
localization error, in a sense that will be described next.
Additional parameters for the optimization problem, largely motivated by our experimental setup, include sphere diameter a = 1.9 cm, vibration axis angle φ = 0◦ , vibration
frequency 40 Hz, and vibration amplitude A = 0.191 cm. In order to ensure good lo49

calization performance when the dipole is at diﬀerent locations within the working area,
we divide the area into a grid of 1 × 1 cm2 cells, resulting in a total of 210 grid points
(x = 0, ±1, · · · , ±10, y = 1, 2, · · · , 10). For a given intra-sensor spacing, for each grid point
(x(k), y(k)), 1 ≤ k ≤ 210, we evaluate the corresponding CRBs Cx (k) and Cy (k) for the
2
2
estimation variances σx (k) and σy (k), assuming that the dipole is located at (x(k), y(k)).

Note that Cx (k) and Cy (k) represent the minimum uncertainty in estimating x(k) and y(k),
respectively. There are two potential ways to deﬁne a single quantity characterizing the
localization uncertainty. The ﬁrst criterion is based on the so-called D-optimality [30] and it
represents the area of the ellipsoidal region of maximum conﬁdence and can be evaluated using the determinant of the corresponding 2 × 2 sub-matrix of I −1 (θ). This criterion could be
misleading, especially when the ellipsoidal region is thin. Therefore, we have adopted an alternative based on A-optimality and capture the bound on the localization uncertainty using
U (k) = π Cx (k) + Cy (k) . The cumulative uncertainty for a given intra-sensor spacing is
then obtained as Uc =

k U (k). Fig. 3.2 shows the cumulative uncertainty Uc as a function

2
of the intra-sensor spacing, for three diﬀerent noise variances σm for the measurement noise.

Here m1 = 0.0065 cm2 /s2 represents a typical value of variance for our IPMC prototype ﬂow
sensors. It is clear that there is an optimal range (around 2 cm) for the intra-sensor distance
where the cumulative localization uncertainty is minimized. Larger or smaller spacing leads
to bigger uncertainty. This can be explained as follows. When the sensors are too packed,
they cannot pick up suﬃciently distinct signals (think about the extreme case of all sensors
stacked together) and cannot cover a wide area with large signal-to-noise ratio. On the other
hand, when the sensors are very far apart, some sensors become too distant from the dipole
source and their signal-to-noise ratios become very low, negatively impacting the localization

50

performance. Another interesting observation from Fig. 3.2 is that, the optimal intra-sensor
spacing is not sensitive to the value of the noise variance. This is a positive news for the
2
design, since it is often diﬃcult to know the precise value of σm in practice.

Uc (cm2)

1.5

σ2 = 0.5 m
m
2
σm
σ2
m

1

1

= m1
= 1.5 m

1

0.5

0
0

1

2

3

4

Intra−sensor separation (cm)

Figure 3.2: CRB-based cumulative uncertainty in localization as a function of intra-sensor
spacing.

3.3.2

Localization of Dipole Source with Unknown Vibration Amplitude and Orientation

We now present simulation results on the performance of the proposed Gauss-Newton and
Newton-Raphson schemes on localizing a dipole source with unknown vibration amplitude
and orientation. Fig. 3.3 illustrates the simulation setup. The lateral line system, consisting
of 6 ﬂow sensors with intra-sensor spacing of 2 cm, is placed parallel to the x−axis and
centered at the origin. The dipole source is placed at 19 diﬀerent locations along an ellipsoidal
51

track that is centered at (0, 6) cm. Speciﬁcally, the dipole locations are prescribed by

(xs (k), ys (k)) = (10 cos ψ(k), 6 + 4 sin ψ(k)) , k = 1, 2, · · · , 19,
where ψ(k) =

(3.24)

(k−1)π
9 . Note that point 1 and point 19 overlap, as shown in Fig. 3.3. In order

to test the algorithms’ capability in estimating vibration amplitude A and orientation φ in
addition to the source location, we vary A (in cm) and φ along the track as follows:

A(k) = 0.191 − 0.01(k − 1),

(3.25)

(k − 1)π
.
9

(3.26)

φ(k) =

2
The sphere diameter a, vibration frequency, and the measurement noise variance σm are set

to be 1.9 cm, 40 Hz, and m1 (0.0065 cm2 /s2 ), respectively. The stopping criterion ǫ for both
the Gauss-Newton algorithm and the Newton-Raphson algorithm is chosen to be 0.01.
Fig. 3.4(a) shows the actual and estimated source locations under the Gauss Newton
(GN) method and the Newton Raphson (NR) method for a typical run. It can be seen the
localization is well accomplished under both schemes. Furthermore, Fig. 3.4(b) shows the
localization error for each of the 19 points on the track, where the error is deﬁned as the
Euclidean distance between the actual and estimated locations. The two proposed methods
have comparable performance in terms of the localization error. The maximum error under
the GN method is 0.38 cm and that under the NR method is 0.39 cm.
Fig. 3.5 shows the estimated source vibration amplitudes and the estimation errors for
both methods. The maximum amplitude error is 0.0013 cm and 0.0015 cm for the GN and
NR methods, respectively. Fig. 3.6 shows the estimated source vibration orientations and
52

12

5

10

y (cm)

8
6

1(19)

10

4
2
0
−15

15
Sensors

−10

−5

0
x (cm)

5

10

15

Figure 3.3: Illustration of the simulation setup, where 19 points as indexed in the ﬁgure will
be localized.
the corresponding errors. The maximum orientation error is 0.12 rad and 0.13 rad for the
GN and NR methods, respectively.
From Figs. 3.4, 3.5, and 3.6, we can see that both the GN and NR methods are able
to simultaneously estimate the location, vibration amplitude, and vibration orientation successfully. In addition, while the GN method has slightly higher estimation accuracy for all
variables, the overall performances of the two methods are quite close.

3.3.3

Further CRB Analysis on Performance

We have further used the CRB analysis to evaluate the localization performance of our
proposed methods. In this evaluation, a similar setup as in Subsection 3.3.2 is considered,
53

Actual locations
Gauss−Newton method
Newton−Raphson method

12
10

y (cm)

8
6
4
2
0

Sensors
−10

−5

0

5

10

x (cm)

Localization error (cm)

(a)

Gauss−Newton method
Newton−Raphson method

0.5
0.4
0.3
0.2
0.1
0
0

5

10

15

20

Index of points
(b)
Figure 3.4: Simulation results: (a) Estimation of source locations; (b) localization errors at
the 19 points.

54

Actual amplitude
Gauss−Newton method
Newton−Raphson method

Amplitude (cm)

0.2
0.15
0.1
0.05

Amplitude error (cm)

0
0

−3

5

10

15

20

5

10

15

20

x 10
1.5
1
0.5
0
0

Index of points
Figure 3.5: Simulation results: Estimation of the source vibration amplitudes.

except that we ﬁx the vibration amplitude and orientation at all points to be 0.191 cm
and 0◦ , respectively, and treat them as known. We run the simulation of localizing the
19 points using the GN and NR algorithms for 20 times, based on which the empirical
2
2
estimation variances σx and σy for xs and ys are obtained. Fig. 3.7 shows the comparison

of the square-roots of empirical variances with the corresponding CRBs. It is clear that the
CRBs are lowest when the source location is closest to the lateral line (points # 13–16),
which is explained the highest signal-to-noise ratio (SNR) at those locations. The empirical
55

Source direction (rad.)
Direction error (rad.)

Actual direction
Gauss−Newton method
Newton−Raphson method
6
4
2
0
0

5

10

15

20

5

10
Index of points

15

20

0.2

0.1

0
0

Figure 3.6: Simulation results: Estimation of the source vibration orientations.

estimation variances under both schemes are higher than the CRBs; however, the diﬀerences
are not big, especially at the points where the SNR is large. This analysis conﬁrms that the
proposed GN and NR schemes deliver close-to-optimal estimation performance.
56

Cramer−Rao bound
Gauss−Newton method
Newton−Raphson method
0.2
0.15

x

σ (cm)

0.25

0.1
0.05
0
0

5

10

15

20

Index of points
(a)

Cramer−Rao bound
Gauss−Newton method
Newton−Raphson method

0.25

σy (cm)

0.2
0.15
0.1
0.05
0
0

5

10

15

20

Index of points
(b)
Figure 3.7: Simulation results: Comparison of square-root of empirical estimation variances
(a) σx and (b) σy with the corresponding CRBs.

57

3.4
3.4.1

Experimental Setup and Flow Model Veriﬁcation
Experimental Setup

We have experimentally examined the performance of the proposed GN and NR algorithms
and the beamforming algorithm reported in the literature, using an artiﬁcial lateral line
prototype consisting of six ionic polymer-metal composite (IPMC) sensors, shown in Fig. 2.4.
An IPMC has three layers, with an ion-exchange polymer membrane sandwiched by metal
electrodes. Inside the polymer, (negatively charged) anions covalently ﬁxed to polymer
chains are balanced by mobile, (positively charged) cations. Deformation under a mechanical
perturbation redistributes the cations, producing a detectable electric signal (e.g., shortcircuit current) that is well correlated with the mechanical stimulus, which explains the
sensing principle of an IPMC [81]. Each sensor in our lateral line prototype, measuring 8 mm
long, 2 mm wide, and 200 µm thick, was cut from an IPMC sheet fabricated by the Smart
Microsystems Laboratory at Michigan State University, following a recipe similar to the one
described in [49]. Based on the optimization results from CRB analysis (Subsection 3.2.4),
the intra-sensor separation was set to be 2 cm, resulting in a total span of 10 cm for the
sensor. The latter can be considered as the Body Length (BL) of the lateral line system,
and we will test its capability in localizing a dipole source within a distance of about one
BL, which is the typical working range for biological lateral lines. All sensors can bend in
the direction that is parallel to the lateral line, thus measuring the ﬂow component in the
same direction.
Fig. 3.8(a) illustrates the experimental setup. Experiments are conducted in a tank
that is 6 feet long, 2 feet wide, and 2 feet deep. The short-circuit current output of each
58

IPMC sensor is ampliﬁed with a two-stage ampliﬁcation circuit, the details of which can
be found in [35]. Acquisition and processing of the IPMC sensor output are conducted
through a dSPACE system (DS1104, dSPACE Inc., Germany). A digital low-pass ﬁlter is
further implemented to remove high-frequency noises from the sensor signals. The dipole
source is created with a mini-shaker (Model 4810, Br¨el & Kjaer, Denmark)(Fig. 2.6(b)), the
u
vibration amplitude and frequency of which can be readily controlled through a voltage input
to the mini-shaker. A lightweight bar ﬁrmly attached to the mini-shaker then translates the
vibration to a sphere rigidly coupled to the bar. The sphere, which is a aluminum ball, has
a diameter of 1.9 cm. The sphere and the IPMC sensors are submerged underwater, about
5 cm from the water surface. The dipole source location and vibration direction with respect
to the IPMC lateral line can be adjusted by moving the stand holding the IPMC lateral line.
The dipole frequency is chosen to be 40 Hz, which is in the typical range for the study of
biological and artiﬁcial lateral lines. The vibration amplitude of the sphere is measured with
a laser displacement sensor (OADM 20I6441/S14F, Baumer Electric).

3.4.2

Flow Model Validation

The basis of our estimation algorithms is the analytical model (3.5) for a dipole sourcegenerated ﬂow ﬁeld. It is thus of interest to examine the validity of this model experimentally.
The validation has been conducted using a digital particle image velocimetry (DPIV) system
(LaVision, Ypsilanti, MI). In a DPIV system, small particles are dispersed in a ﬂuid and a
laser sheet is created in the ﬂuid to illuminate the particles. Processing of images taken in
quick successions can reveal the movement of particles and thus provide information about
the ﬂow ﬁeld.
59

Power
supply
dSPACE
Vibrating sphere
dipole
Electronic
conditioning

Minishaker
Supporting
frame

IPMC
Sensors

IPMC bending
direction

Tank

Figure 3.8: Experimental setup: Schematic of the experimental system.

We measured the x−component of the ﬂow velocity at the origin when a dipole source,
enabled with a slightly diﬀerent mechanism from that in Fig. 3.8(b) [3], was placed at
eight diﬀerent locations: (±1.5, 1.5) cm, (±3, 1.5) cm, (±1.5, 3) cm, and (±3, 3) cm. Since
the DPIV system has a measurement frequency of 15 Hz, we have obtained the velocity
amplitude by taking the maximum measurement over a suﬃciently long period (over 10 s).
Fig. 3.9 shows the comparison of the DPIV measurement with the prediction from the
theoretical model (3.5). As can be seen, the match between the model and the experimental
60

y = 1.5 cm: model
y = 1.5 cm: measurement
y = 3 cm: model
y = 3 cm: measurement

Flow velocity (cm/s)

6
5
4
3
2
1
0
−4

−2

0

2

4

x (cm)
Figure 3.9: Experimental results: Comparison of DPIV-measured ﬂow velocity with the
model prediction.
measurement is good, with the maximum error being 0.2 cm/s among the eight sampled
points.

3.4.3

Calibration of IPMC Sensors

The relationship between the current output I of an IPMC sensor and the local ﬂow velocity
v can be captured by the cascaded connection of ﬂuid-structure interaction dynamics G1 and
IPMC transduction dynamics G2 , both of which can be treated as linear when the ﬂow is
relatively small. The dynamics G1 describes how the tip displacement q of the sensor beam
responds to the ﬂow v, and can be modeled following the approach presented by McHenry et
al. in modeling the interactions between a superﬁcial neuromast and the ambient ﬂow [61].
61

The dynamics G2 , on the other hand, relates the current output I to the tip displacement
q and can be captured with a physics-based model [16]. It can be shown that, based on the
physical parameters for the IPMC sensors and the operating frequency in our experiments,
the transfer function G1 (s) can be approximated by [61]

G1 (s) =

c
q(s)
= 1,
v(s)
s

(3.27)

and G2 by [13]

G2 (s) =

I(s)
= c2 s,
q(s)

(3.28)

where s is the Laplace variable, and c1 and c2 are lumped physical constants. The detailed
derivation and approximation of G1 and G2 are explained earlier in Chapter 2.

From (3.27) and (3.28), the short-circuit current output of an IPMC sensor is (approximately) proportional to the instantaneous local ﬂow velocity. The proportionality constant
can vary from sensor to sensor because of the spatial inhomogeneity in the fabricated IPMC
material, diﬀerence in exact sensor dimensions, and discrepancy across diﬀerent channels
of the sensing circuit. The constant is identiﬁed through a calibration procedure described
next. We put the mini-shaker-based dipole source at diﬀerent locations with respect to each
sensor and extracted the amplitude of sensing current with FFT. The theoretical value of
the ﬂow velocity at the location of the sensor was obtained using (3.5). Fig. 3.10 shows
the IPMC signal amplitude versus the theoretical local velocity for one IPMC sensor, from
which an approximately linear relationship between the sensor output I and the ﬂow velocity
I
v parallel to the lateral line can be observed. The proportional constant η = v is identiﬁed

62

using the Matlab command polyf it(·, ·, 1).
The proportionality constants identiﬁed for the six IPMC sensors are η1 = 0.49 µA/(cm · s−1 ),
η2 = 0.21 µA/(cm · s−1 ), η3 = 0.165 µA/(cm · s−1 ), η4 = 0.23 µA/(cm · s−1 ), η5 =
0.24 µA/(cm · s−1 ), and η6 = 0.3 µA/(cm · s−1 ). Given ηi , the ﬂow velocity at the site
of sensor i is calculated from the sensor output Ii via

IPMC current (µ A)

I
vi = i .
ηi

Measurement
Linear fit

0.25
0.2
0.15
0.1
0.05
0
0

0.1

0.2

0.3

0.4

0.5

Flow velocity (cm/s)
Figure 3.10: Illustration of IPMC sensor calibration.

63

0.6

3.4.4

Experimental Results on Dipole Source Localization

We have conducted localization experiments in a very similar setup as for the simulation (see
Fig. 3.3). Speciﬁcally, we have placed the dipole source (relative to the lateral line) at 19
points along an ellipsoidal track centered at (0, 6) cm, as prescribed by (3.24). The vibration
amplitude at each point varies according to (3.25). Due to the diﬃculty in precisely orienting
the dipole vibration axis, we have chosen to use only two diﬀerent orientations, 0 and π rad.
2
In particular, we set
φ(k) =



 0, k = 1, · · · , 10

 π
 , k = 11, · · · , 19
2

For each location, sensor outputs of 10 s are acquired at a sampling rate of 1 kHz. FFT
is performed on the data to extract the signal amplitudes, which are then used for estimation with the GN and NR methods. For comparison purposes, the same (raw) sensor
data have been used for estimation using the beamforming algorithm; see the description in
Subsection 3.2.3.
Fig. 3.11 shows the experimental results on the localization. It is notable that the localization performance has only slightly degraded from that in simulation (Fig. 3.4). The
maximum localization errors for the GN and NR methods are 0.46 cm and 0.48 cm, respectively, which occur at point #19, where the vibration amplitude is the smallest. Fig. 3.12(a)
and (b) show the estimated source vibration amplitudes and source vibration directions at
the 19 points. For both methods, the maximum errors in the amplitude estimation are
around 0.002 cm, about 50% larger than those in simulation. The maximum errors in the
orientation estimation are around 0.12 rad, which are almost the same as those in simulation;
this might be explained by the fact that only two diﬀerent orientations were tested in the
64

experiments. Another observation is that, in consistency with the simulation results, the
GN method has demonstrated slightly better performance than the NR method.

Actual locations
Gauss Newton method
Newton−Raphson method

12
10

y (cm)

8
6
4
2
0

−10

−5

Sensor
0

5

10

x (cm)

Localization error (cm)

(a)

Gauss−Newton
Newton−Raphson

0.5
0.4
0.3
0.2
0.1
0
0

5

10

15

20

Index of points
(b)
Figure 3.11: Experimental results: (a) Estimation of source locations; (b) localization errors
at the 19 points.
We have further examined the performance of the beamforming algorithm on the same
65

0.2

0.1

Amplitude error (cm)

Amplitude (cm)

Actual amplitude
Gauss−Newton method
Newton−Raphson method

0
0 −3
x 10
3

5

10

15

20

5

10
Index of points

15

20

2
1
0
0

Direction error (rad.) Source direction (rad.)

(a)
Actual direction
Gauss−Newton
Newton−Raphson

2

1

0
0
0.2

5

10

15

20

5

10
Index of points

15

20

0.1

0
0

(b)
Figure 3.12: Experimental results: (a) Estimation of the source vibration amplitudes; (b)
estimation of the source vibration orientations.

66

data. In order to create the energy-like map for each dipole state, one needs to discretize the
space for the dipole state (location, vibration amplitude, and vibration orientation). Clearly,
the resolution of the discretization has a direct impact on the estimation performance; a
ﬁner grid is expected to produce more accurate resolution (sharper map), but at the cost of
computational time. To ease the computation, we have assumed that the vibration amplitude
and orientation at each location are known, and thus focused on only the estimation of
source location. We have used three diﬀerent resolutions for discretizing the working area,
0.02 × 0.02 cm2 , 0.05 × 0.05 cm2 , and 0.1 × 0.1 cm2 . The dipole location is estimated by
visually ﬁnding the maximum of the energy-like map. As an illustration, Fig. 3.13 shows
the energy-like map when the dipole source was placed at point #3, with the 0.1 × 0.1 cm2
discretization grid, where the peak is clearly visible.

60

E

40
20
0
100

100
0

50

y (mm)

x (mm)
0 −100

Figure 3.13: Experimental results: Constructed energy-like map when the dipole is located
at point #3 (0.1 × 0.1 cm2 grid).
Fig. 3.14 shows the localization results corresponding to the three discretization schemes
67

Table 3.1: Experimental results: maximum and average estimation errors under diﬀerent
algorithms. “BF (0.02)” stands for beamforming with 0.02 × 0.02 cm2 grid, and similar
interpretations apply for “BF (0.05)” and “BF (0.1)” . Note that vibration amplitude and
orientation are assumed known for the BF algorithms.
Algorithms

GN

NR

BF (0.02)

BF (0.05)

BF (0.1)

Average localization error (cm)
Max. localization error (cm)
Average amplitude error (cm)
Max. amplitude error (cm)
Average orientation error (rad)
Max. orientation error (rad)

0.21
0.46
0.0014
0.0019
0.1
0.12

0.21
0.48
0.0014
0.0021
0.1
0.14

0.48
0.91
-

0.65
1.1
-

0.89
1.33
-

under the beamforming approach. As one would expect, the localization error decreases as
the grid gets ﬁner; however, the error does not scale linearly with the grid resolution. The
maximum localization errors for the three discretization schemes (ﬁner to coarser) are 0.91
cm, 1.1 cm, and 1.33 cm, respectively. For a more clear comparison, Table 3.1 summarizes
the average and maximum localization errors for the 19 points, for the Gauss-Newton (GN),
Newton-Raphson (NR), and three beamforming (BF) algorithms. It can be seen that, with
the ﬁnest grid (0.02 × 0.02 cm2 ), the BF algorithm results in an average/maximum localization error that is twice as much as what is resulted from the GN and NR algorithms.

To put the comparison in perspective, we have further recorded the computation time
taken by each algorithm. All algorithms have been run on Samsung laptop with a 2.53 GHz
processor and 4 GB RAM. The localization problem at all 19 points is executed under each
algorithm for seven times, and Table 3.2 summarizes the average time per point taken by
each algorithm. Note that for GN and NR, the computation includes FFT and executing
the recursive update until convergence, and for the beamforming algorithm, the computation
68

Actual track
BF (0.02 × 0.02)
BF (0.05 × 0.05)
BF (0.1 × 0.1)
10

y (cm)

8
6
4
2
0

−10

−5

Sensors
0

5

10

x (cm)

Localization error (cm)

(a)

0.02 × 0.02 lattices
0.05 × 0.05 lattices
0.1 × 0.1 lattices

1.5

1

0.5

0
0

5

10

15

20

Index of points
(b)
Figure 3.14: Experimental results: Localization of the dipole source using the beamforming
algorithm: (a) Estimation of the source locations; (b) localization errors at the 19 points.

69

Table 3.2: Average computational time per localization point for diﬀerent algorithms.
Algorithms

Average time (s)

GN
NR
BF (0.02 × 0.02 cm2 )
BF (0.05 × 0.05 cm2 )
BF (0.1 × 0.1 cm2 )

1.1
1.2
11.2
2.1
1.26

includes evaluating the correlation matrix and constructing the energy-like map. From Table 3.2, it can be seen that GN is slightly more eﬃcient than NR. With a discretization grid
of 0.02 × 0.02 cm2 , the beamforming algorithm takes a comparable amount of computation
time as GN or NR; however, from Table 3.1, the corresponding maximum localization error
is three times as much. Considering that the GN and NR algorithms simultaneously estimate four variables while the beamforming algorithm only does two (in this setup), we can
conclude that the proposed GN and NR algorithms are more accurate and computationally
eﬃcient than the beamforming algorithm. One likely reason that the beamforming algorithm
produces larger estimation error is the noise in the raw signals from IPMC sensors, which
would have a lesser impact for the GN and NR schemes because of the FFT procedure.

3.5

Chapter Summary

In this chapter, using the analytical model for a dipole-generated ﬂow ﬁeld, we formulate
a nonlinear estimation problem and present two novel iterative schemes for simultaneously
localizing a dipole source and estimating its vibration amplitude and orientation. The ﬁrst
scheme, which is based on the Gauss-Newton (GN) method, solves the nonlinear estimation
problem through iterative linearization. In the second scheme, the Newton-Raphson (NR)
70

method is used to solve the nonlinear equation resulting from the ﬁrst-order optimality
condition. The ﬂow model, with measurement noise properly incorporated, is also used to
derive the Cramer-Rao bound (CRB). Analysis based on the CRB is subsequently exploited
for the optimal design of intra-sensor spacing of the lateral line. Simulation is conducted
to localize a dipole source placed at diﬀerent points along an ellipsoidal track, with its
vibration amplitude and orientation varying from one point to next. Simulation results have
shown that both the GN and NR algorithms are able to simultaneously estimate the source
location, vibration amplitude and orientation with comparable precision. Furthermore, the
comparison with the CRB shows that the algorithms achieve near-optimal performance at
many points.
We have further validated the algorithms using an experimental prototype of lateral line
that is comprised of six millimeter-scale ionic polymer-metal composite (IPMC) ﬂow sensors.
With excellent agreement with the simulation results, the experimental results conﬁrm that
the GN and NR schemes have comparable performance in accuracy and computational eﬃciency with the GN scheme having a slight advantage. Speciﬁcally, for the body length (BL)
of 10 cm, the maximum localization error is less than 5% of BL when the source is within
the distance of one BL. Experimental results have also shown that the proposed schemes
are superior to the beamforming method; when the beamforming method is only required
to determine the dipole location (with vibration amplitude/orientation given), its maximum
(average, resp.) localization error is three (four, resp.) times of those under the GN and NR
schemes while consuming comparable computing time.

71

Chapter 4
Localization of a Moving Dipole
Source
Existing studies on artiﬁcial lateral lines (ALLs) have been mostly focused on the localization
of a ﬁxed underwater vibrating sphere (dipole source). In this chapter we introduce a novel
analytical model for the problem of localizing a moving dipole source using an ALL. Such
a moving source could represent a swimming ﬁsh or robotic ﬁsh, which demonstrates a
combination of translational motion and oscillatory body/ﬁn motion.
First, we formulate a nonlinear estimation problem based on an analytical model for
the moving dipole-generated ﬂow ﬁeld, where we assume that the source location, vibration
amplitude, and moving speed are unknown. We also assume that the amplitudes of the ﬂow
velocities at the DC and dipole vibration frequencies are available. Such amplitudes can
be obtained by conducting sliding digital Fourier transform (SDFT) of the measured signal
in the time domain. The estimation problem is solved recursively with the Gauss-Newton
method, and we illustrate the proposed approach with simulation and experimental results.
72

The remainder of the chapter is organized as follows. The localization of the traveling
dipole source is formulated in Section 4.1. The sliding discrete Fourier transformation algorithm is described in Section 4.2, and the Gauss-Newton recursive estimation algorithm is
presented in Section 4.3. Simulation and experimental results are presented in Section 4.4
and Section 4.5, respectively. Finally, concluding remarks are provided in Section 4.6.

4.1

Localization of a Traveling Dipole Source: Problem
Formulation

We now formulate the problem of localizing a moving dipole source in the two-dimensional
x − y plane. Consider Fig. 4.1, where an ALL comprising N sensors is located in parallel
to the x-axis, with the sensor locations denoted as (xi , yi ), 0 ≤ i ≤ N − 1. We denote the
location of a dipole source at time t by (xs (t), ys (t)). The source is assumed to move with a
constant speed vm parallel to the ALL, and vibrating with a velocity v d with the vibration
axis perpendicular to the ALL. We further assume that the presence of sensors has negligible
eﬀect on the ﬂow distribution generated by the moving source. Each sensor is assumed to
provide a (noisy) measurement of the local ﬂow velocity v along the x−direction.
Under an assumption of ideal ﬂows, the ﬂow velocity (x-component) at the site of sensor
i can be derived as [53]

v(xi , yi , t) =

a3
(2(xi − xs (t))2 − (yi − ys (t))2 vm
5
2 ri (t)
+3(xi − xs (t))(yi − ys (t))v d (t) ,

(4.1)

where ri (t) = (xs (t) − xi , ys (t) − yi )T , a represents the radius of the sphere, and T denotes
73

Figure 4.1: Illustration of the problem setup for the localization of a moving dipole source.

the transpose. Assume that the sphere vibrates with angular frequency ω and amplitude A:

v d (t) = A sin(ωt).

When the source location (xs (t), ys (t)) is changing relatively slowly with respect to the dipole
motion v d (t), Eq. (4.1) can be viewed as the sum of a constant term and an oscillatory term,
corresponding to the source translation and vibration motions, respectively. Denote the
amplitudes of the two terms as m1,i and m2,i , respectively; namely,
a3
(2(xi − xs )2 − (yi − ys )2 )vm ,
5
2 ri
a3
△
(3(xi − xs )(yi − ys )A),
m2,i =
2 ri 5
△

m1,i =

(4.2)
(4.3)

where θ = (vm , A, xs , ys ) represents the set of parameters to be estimated. Note that θ
will have a dependence on t since the source location (xs (t), ys (t)) changes over time, the
74

magnitude of these components namely f1,i = m1,i , and f2,i = m2,i . We assume that
noise-corrupted measurements of m1,i , m2,i , i = 0, 1, · · · , N − 1, are available:
m1,i =
ˆ

m1,i + d1,i ,

m2,i =
ˆ

m2,i + d2,i ,

N −1
where d1,i and d2,i are noises associated with the sensor i. We assume the noises {d1,i , d2,i }i=0

are of zero mean and not correlated with each other. We further deﬁne

△

m1 =
△

m2 =
△

m1 =
ˆ
△

m2 =
ˆ
△

f1 (θ) =
△

f2 (θ) =

m1,0 , · · · , m1,N −1

T

,

m2,0 , · · · , m2,N −1

T

,

m1,0 , · · · , m1,N −1
ˆ
ˆ

T

,

m2,0 , · · · , m2,N −1
ˆ
ˆ

T

,

f1,0 (θ), · · · , f1,N −1 (θ)

T

,

f2,0 (θ), · · · , f2,N −1 (θ)

T

,

The estimation problem is formulated as follows. Let the sphere diameter a, the sensor
N −1
locations {(xi , yi )}i=0 , and the frequency ω be known. Given the sensor measurements m1
ˆ

ˆ
and m2 at the time t, provide an estimate of the parameter θ(t), θ(t) = (ˆm , A, xs (t), ys (t)),
ˆ
v ˆ ˆ
ˆ
such that

ˆ
ˆ
ˆ
ˆ
ˆ
J(θ) = β(m1 − f1 (θ))T (m1 − f1 (θ)) + (1 − β)(m2 − f2 (θ))T (m2 − f2 (θ))
ˆ
ˆ
ˆ
ˆ

75

(4.4)

is minimized, where β ∈ [0, 1] is a parameter weighing the importance of the DC term relative
to the oscillatory term.

4.2

The Sliding Discrete Fourier Transform (SDFT)
Algorithm

Since the actual sensor output is a discrete-time signal, discrete Fourier transform (DFT) can
be used to extract the signal amplitudes at diﬀerent frequencies. However, because the target
is moving, the amplitudes m1 and m2 evolve with the time index n and it is thus not practical
to collect a large sequence of sensor data and perform DFT to extract the measurements
m1 and m2 . The SDFT algorithm [44, 82] is particularly suitable for computing a speciﬁc
ˆ
ˆ
spectral bin in real-time application based on a sliding window of time samples. For this
reason, the SDFT algorithm has been adopted to compute the two frequency components of
interest in this paper.
For a window of M time samples, {x [n − (M − 1)] , · · · , x [n]}, we can evaluate the p−th
spectral bin using DFT:

Xp [n] =

M −1
i=0

x [n − M + 1 + i] e−2jπpi/M , for 0 ≤ p ≤ M − 1.

(4.5)

Recall (4.1). In our target localization problem, there are two frequencies that are of interest,
the DC signal (spectral bin p = 0) and the signal with angular frequency ω (spectral bin
p =

ωM
2πFs ), where [·] denotes rounding to the nearest integer and Fs is the sampling

frequency for obtaining sensors signals x [·]. Exploiting the fact that all elements but one
in neighboring windows are identical, the SDFT algorithm can eﬃciently evaluate Xp [n]
76

recursively:
Xp [n] = e−jω0 (Xp [n − 1] − x [n − N ] + x [n]),

(4.6)

where ωo = 2πp/M .

For p = 0, from (4.5), we can see that X0 [n] essentially sums up the M signal samples
in the window. Consequently, the magnitude of the DC component for {x[·]} is

|X0 [n]|
M . For

ωM
the AC component (p = 2πF ) is given as
s

Xp [n] =

M −1
i=0

m2 sin(wo (n − M + 1 + i))e−2jπpi/M

m e−jwo (M −n)
= 2
2j
=
=

M −1
i=0

m2 e−jwo (M −n)

(1 − e−2jwo i )

M−

2j
m2 e−jwo (N −n)
M
2j

1 − e−2jwo M
1 − e−2jwo

(4.7)

where e−2jwo M − 1 = 0, because p is an integer number. Taking the norm for both sides:
Xp [n] =

m2 M
2

(4.8)

2Xp [n]
M

(4.9)

Then, the corresponding AC component is

f2 = |m2 | =
77

4.3

The Gauss-Newton Method for Solving the Estimation Problem

While f1 (θ) and f2 (θ) are linear in vm and A, they are highly nonlinear in (xs , ys ). To
eﬃciently solve the estimation problem, we adopt the Gauss-Newton method. First we
¯
approximate f1 and f2 by linearizing them around some nominal point θ:

¯
¯
¯
¯
f1 (θ) = f1 (θ) + B1 (θ)(θ − θ),
¯
¯
¯
¯
f2 (θ) = f2 (θ) + B2 (θ)(θ − θ),
where
∂f1
∂f2
¯
¯
B1 (θ) =
, B2 (θ) =
.
∂θ θ=θ
∂θ θ=θ
¯
¯
Accordingly, we modify the cost function for the estimation problem from (4.4) to

ˆ
¯ ˆ
¯ ˆ
¯ ˆ
¯ ˆ
J(θ) = β(m1 − f1 (θ))T (m1 − f1 (θ)) + (1 − β)(m2 − f2 (θ))T (m2 − f2 (θ)),
ˆ
ˆ
ˆ
ˆ

(4.10)

which becomes a standard least squares estimation problem. A solution to this problem is

ˆ
¯
¯
¯
¯
¯
θ = θ + λ1 βB1 (θ)T B1 (θ) + (1 − β)B2 (θ)T B2 (θ)

−1

T ¯ ˆ
T ¯ ˆ
¯
¯
βB1 (θ)(m1 − f1 (θ)) + (1 − β)B2 (θ)(m2 − f2 (θ)) ,

(4.11)

where λ1 is a stepping parameter for the DC component.
In order to minimize the impact of error introduced by linearization, we use a recursive
78

version of (4.11) to update the estimate:

ˆ
θk+1 =

ˆ
ˆ
ˆ
ˆ
ˆ
θk + λ1 βB1 (θk )T B1 (θk ) + (1 − β)B2 (θk )T B2 (θk )

−1

T ˆ
T ˆ ˆ
ˆ
ˆ
βB1 (θk )(m1 − f1 (θk )) + (1 − β)B2 (θk )(m2 − f2 (θk )) ,
ˆ

(4.12)

ˆ
where the initial estimate θ0 is chosen properly. When β = 0 (i.e., considering the oscillatory
component only), the algorithm (4.12) is simpliﬁed to

ˆ
ˆ
ˆ
ˆ
θk+1 = θk + λ1 B2 (θk )T B2 (θk ))

−1

T ˆ
ˆ
ˆ
B2 (θk )(m2 − f2 (θk )) .

(4.13)

ˆ
ˆ
The iteration is stopped when θk+1 − θk ≤ ǫ, where ǫ > 0 is a speciﬁed tolerance.

4.4

Simulation Results

Fig. 4.2 illustrates the simulation setup. The lateral line system is placed parallel to the
x−axis and centered at (0, 5) cm. It consists of 6 sensors, with the sensor-to-sensor separation
of 2 cm. The dipole source (sphere diameter 1.9 cm) vibrates at 40 Hz with amplitude
A = 0.191 cm. The source moves from left to right with a constant speed of 1.5 cm/s. The
initial location of the source is (0, 2) cm, and the time duration for the movement is 6.667 s,
resulting in a terminal location of (10, 2) cm.
Fig. 4.3 shows the simulated ﬂow velocities (x−component) at each sensor site. Figs. 4.4
and 4.5 show the magnitudes of the DC and AC components, respectively, obtained through
SDFT. The sampling frequency used is 20 kHz, and the sliding window size is 500 samples
(25 ms).
79

(x (t ), y (t ))
1

1

(x (t ), y (t ))
n

(x (t ), y (t ))

n

f

f

2

y (cm)

vm=1.5 cm/s

Sensors
0

2

4

6

8

10

12

x (cm)

Figure 4.2: Illustration of the simulation setup.

It is clear from Figs. 4.4 and 4.5 that both the DC and AC components contain relevant
information for localizing the source, as can be seen from the shifting magnitude proﬁles in
the ﬁgures. We assume a zero-mean Gaussian noise with a variance of 0.0065 (cm2 /s2 ) for
the sensors.
It is of interest to investigate how the weighting parameter β in (4.4) impacts the estimation performance. Fig. 4.6 shows the actual and estimated source locations and the tracking
error at each point under the proposed method for a typical run, where the stopping criterion
ǫ is chosen to be 0.01, and at period time ∆ = 0.2 s over the period (10 s), we have estimated
the source locations and the other parameters. The maximum tracking errors are 0.044 cm
for β = 0 and 0.0441 for β = 0.5, respectively. Table 4.1 summarizes the averages of the
estimates for other parameters (vm and A), where each average is obtained by taking the
mean of the estimated values over the period [0.2, 10] s. When β = 0, we do not have direct
80

(a)

(b)

10
0
−10
−20
0

20

Sensor 3(cm/s)

20

Sensor 2(cm/s)

Sensor 1(cm/s)

20

10
0
−10
−20
0

5

−20
0

5

5

20

Sensor 6(cm/s)

Sensor 5(cm/s)

−10

−10

Time (s)
(f)

20

0

0

Time (s)
(e)

20
10

10

−20
0

5

Time (s)
(d)
Sensor 4(cm/s)

(c)

10
0
−10
−20
0

5

Time (s)

Time (s)

10
0
−10
−20
0

5

Time (s)

Figure 4.3: Simulation results: Time trajectories of the ﬂow velocities at the sensor locations.

information to estimate the moving velocity, so we have used an indirect method to do so
by deﬁning vm as:
¯
x − xk−1
ˆ
ˆ
vm,k = k
¯
,
∆
where k is the discrete time index.

(4.14)

It is clear from Table 4.1 that, the diﬀerence between

vm and vm is small, which will be our justiﬁcation for using just the AC-component in the
¯
ˆ
experimental section.
We have further examined the case where the amplitude of vibration varies while the
81

(a)

(b)

0.6
0.4
0.2
0
0

0.8

Sensor 3(cm/s)

0.8

Sensor 2(cm/s)

Sensor 1(cm/s)

0.8

0.6
0.4
0.2
0
0

5

5

5

Time (s)
(f)
0.8

Sensor 6(cm/s)

Sensor 5(cm/s)

Sensor 4(cm/s)

0.2

0.2
0
0

0.8

0.4

0.4

Time (s)
(e)

0.8
0.6

0.6

5

Time (s)
(d)

0
0

(c)

0.6
0.4
0.2
0
0

0.6
0.4
0.2

5

Time (s)

Time (s)

0
0

5

Time (s)

Figure 4.4: Simulation results: Magnitude trajectories of the DC component at the sensor
sites.

sphere is moving. Fig. 4.7 shows the actual and estimated locations under a typical run,
where the moving velocity is 3 cm/s and the vibration amplitude changes with the time as
shown in Fig. 4.8. From Figs.

4.7 and 4.8, the proposed scheme can eﬀectively estimate

both the location and the vibration amplitude of the dipole. The estimated dipole traveling
speed is vm = 2.87 cm/s, which is close to the actual value of 3 cm/s.
¯
82

(b)

(a)
Sensor 2(cm/s)

10
5

5

Sensor 5(cm/s)

15
10
5
5

Time (s)

15
10
5
0
0

5
Time (s)
(e)

20

20

Sensor 4(cm/s)

10

0
0

5

Time (s)
(d)

0
0

15

5
Time (s)
(f)

20
Sensor 6(cm/s)

Sensor 1(cm/s)

15

20
Sensor 3(cm/s)

20

20

0
0

(c)

15
10
5
0
0

5
Time (s)

15
10
5
0
0

5
Time (s)

Figure 4.5: Simulation results: Magnitude trajectories of the AC component at the sensor
sites.

4.5
4.5.1

Experimental Setup and Results
Experimental Setup

The artiﬁcial lateral line prototype consists of six IPMC sensors as shown in Fig. 4.9. Each
sensor has dimensions 8 mm × 2 mm × 200 µm. The sensor-to-sensor separation was
2 cm, resulting in a total span of 10 cm. Each sensor is mounted in such a way that it would
bend in the direction parallel to the sensor array, which is denoted as x−direction.
83

Actual locations
Estimated locations β = 0
Estimated locations β = 0.5

12

x (cm)

10
8
6
4
2
0
0

1

2

1

2

3

4

5

6

7

3

4

5

6

7

2.5

y (cm)

2
1.5
1
0.5
0
0

Tracking error (cm)

Time (s)
β=0
β = 0.5

0.06
0.04
0.02
0
0

1

2

3

4

5

6

7

Time (s)

Figure 4.6: Simulation results on the estimation of the source location: (a) Localization; (b)
tracking error.

84

Table 4.1: Simulation results: Estimated velocity (vm ), vibrating amplitude (A).
Parameter

Actual

vm (β = 0)
¯

vm (β = 0.5)
¯

vm (β = 0.5)
ˆ

Velocity (cm/s)
Amplitude (mm)

1.5
1.91

1.41
1.87

1.414
-

1.46
1.86

Fig. 4.10 shows the schematic of the experimental setup, while Fig. 4.11 shows the picture of the actual experimental system. The experiments were conducted in a water tank
measuring 6 ft long, 2 ft wide, and 2 ft deep. A custom-built, DC motor-driven conveyor
was used to move a moving dipole in the water. The conveyor track was separated from
the lateral line system by 2 cm. The signals from the IPMC sensors were conditioned and
ampliﬁed before the acquisition and processing by a PC equipped with dSPACE (DS 1104,
dSPACE).

The relationship between the current output I of an IPMC sensor and the

local ﬂow velocity of the AC component m2 can be captured by the cascaded connection
of ﬂuid-structure interaction dynamics G1 and IPMC transduction dynamics G2 , both of
which can be treated as linear when the ﬂow-induced beam deformation is relatively small.
The short-circuit current output of an IPMC sensor is (approximately) proportional to the
instantaneous local ﬂow velocity; the detailed derivation and approximation of G1 and G2
are explained earlier in Chapter 2. The proportionality constant can vary from sensor to
sensor because of the spatial inhomogeneity in the fabricated IPMC material, diﬀerence in
exact sensor dimensions, and discrepancy across diﬀerent channels of the sensing circuit.
The constant is identiﬁed through a calibration procedure described next. We put the minishaker-based dipole source at diﬀerent locations with respect to each sensor and extracted
the amplitude of sensing current with DFT. The theoretical value of the ﬂow velocity at
the location of the sensor is obtained using (4.3). Fig. 4.12 shows the IPMC signal ampli85

tude versus the theoretical local velocity for one IPMC sensor, from which an approximately
linear relationship between the sensor output I and the ﬂow velocity f2 parallel to the
lateral line can be observed. The proportional constant η = fI is identiﬁed using the Mat2
lab command polyf it(·, ·, 1).

The proportionality constants identiﬁed for the six IPMC

sensors are η1 = 0.4 µA/(cm · s−1 ), η2 = 0.25 µA/(cm · s−1 ), η3 = 0.22 µA/(cm · s−1 ),
η4 = 0.33 µA/(cm · s−1 ), η5 = 0.28 µA/(cm · s−1 ), and η6 = 0.23 µA/(cm · s−1 ). Given ηi ,
the ﬂow velocity at the site of sensor i is calculated from the sensor output Ii via
I
vi = i .
ηi
.

4.5.2

Experimental Results on Moving Dipole Source Localization

For a typical run, the dipole source vibrates at 40 Hz with amplitude A = 0.191 cm. The
source moves from left to right with a constant speed of 1.5 cm/s. The initial location of the
source is (0, 2) cm, and the time duration for the movement is 6.667 s, resulting in a terminal
location of (10, 2) cm. During experiments both the IPMC sensors and the vibrating sphere
are completely submerged underwater with a depth of approximately 5 cm. Fig. 4.13 shows
the measured ﬂow velocities (x−component) at each sensor site.
Fig. 4.14 and 4.15 show the experimental magnitudes of the DC and AC components,
respectively, obtained through SDFT. The sampling frequency used is 20 kHz, and the
sliding window size is 500 samples (25 ms).

It is clear from Fig. 4.14 that the magnitude

trajectories of the DC ﬂow component at the sensor sites are subject to uncertain drifts.
Therefore, in the following experimental localization results, we take β = 0 so that only the
86

AC component is used.
Fig. 4.16 shows the actual and estimated source locations and the tracking error at each
point under the Gauss Newton method for a typical run. Table 4.2 lists the estimated moving
speed and vibrating amplitude. The maximum tracking error was 0.33 cm and the errors
in estimating vm and A were less than 10%, which demonstrates the eﬀectiveness of the
proposed estimation approach.
Table 4.2: Experimental results: Estimated velocity (vm ), vibrating amplitude (A).
Parameter

Estimated

Error

Velocity (cm/s)
Vibrating amplitude (mm)

4.6

Actual
1.5
1.91

1.4
2.1

0.1
0.19

Chapter Summary

In this chapter we investigated an analytical model-based estimation approach to the localization of a moving dipole source using an artiﬁcial lateral line consisting of an array of ﬂow
sensors. The moving speed, vibration amplitude, and location of the moving source were
assumed to be unknown, and we proposed a Gauss-Newton algorithm to solve the nonlinear estimation problem based on the evolving signal magnitudes obtained with SDFT. The
eﬀectiveness of the approach was illustrated with simulation and experimental results.

87

Actual locations
Estimated locations β = 0

12

x (cm)

10
8
6
4
2
0
0
2.5

0.5

1

1.5

2

2.5

3

3.5

0.5

1

1.5

2

2.5

3

3.5

2.5

3

3.5

y (cm)

2
1.5
1
0.5
0
0

Time (s)
(a)

Tracking error (cm)

0.03
0.025
0.02
0.015
0.01
0.005
0
0

0.5

1

1.5

2

Time (s)
(b)
Figure 4.7: Simulation results on the estimation of the source location, vibration amplitude
changing with the time: (a) Localization; (b) tracking error.
88

Actual amplitude
Estimated amplitude

0.18
0.16

Amplitude error (cm)

Amplitude (cm)

0.2

6

0
0.5
−3
x 10

1

1.5

2

2.5

3

3.5

1.5
2
Time (s)

2.5

3

3.5

4
2
0
0

0.5

1

Figure 4.8: Simulation results of the source vibration amplitude.

IPMC Sensors

2 cm

Figure 4.9: An experimental prototype of IPMC-based lateral line system.
89

Power
supply
dSPACE
Vibrating sphere
dipole
DC Motor driving
the dipole

Electronic
conditioning

Supporting
frame
IPMC bending
direction

IPMC
Sensors

Moving
direction

Tank

Figure 4.10: Experimental setup of the moving dipole source: Schematic.

90

DC Motor

Mini−shaker
Vibrating
sphere

Sensors

Moving
direction

Figure 4.11: Experimental setup of the moving dipole source: Picture of the experimental
setup.

IPMC current (µ A)

0.25
Measurment
Linear fit

0.2
0.15
0.1
0.05
0
0

0.1

0.2

0.3

0.4

0.5

Flow velocity (cm/s)

Figure 4.12: Calibration of an IPMC sensor.

91

0.6

10
0
−10
0

5

10
0
−10
0

5

−10
0

5

0
−10
0

10
0
−10
0

5

Time (s)

Time (s)

5

Time (s)
(f)
Sensor 6(cm/s)

0

10

Time (s)
(e)
Sensor 5(cm/s)

Sensor 4(cm/s)

Time (s)
(d)
10

(c)
Sensor 3(cm/s)

(b)
Sensor 2(cm/s)

Sensor 1(cm/s)

(a)

10
0
−10
0

5

Time (s)

Figure 4.13: Experimental results: Time trajectories of the ﬂow velocities at the sensor
locations.

92

(a)

(b)

0.5

0
0

1

Sensor 3(cm/s)

1

Sensor 2(cm/s)

Sensor 1 (cm/s)

1

0.5

0
0

5

Time (s)
(e)

5

5

Time (s)
(f)
1

Sensor 6(cm/s)

1

Sensor 5(cm/s)

Sensor 4(cm/s)

1

0.5

0.5

0
0

5

Time (s)
(d)

0
0

(c)

0.5

0
0

5

Time (s)

Time (s)

0.5

0
0

5

Time (s)

Figure 4.14: Experimental results: Magnitude trajectories of the DC component at the
sensor sites.

93

(b)

10
5
0
0

15
10
5
0
0

5

0
0

5

5

15
10
5
0
0

5

Time (s)

Time (s)

5

Time (s)
(f)
Sensor 6(cm/s)

5

10

Time (s)
(e)
Sensor 5(cm/s)

Sensor 4(cm/s)

10

15

0
0

5

Time (s)
(d)
15

(c)
Sensor 3(cm/s)

15

Sensor 2(cm/s)

Sensor 1(cm/s)

(a)

15
10
5
0
0

5

Time (s)

Figure 4.15: Experimental results: Magnitude trajectories of the AC component at the
sensor sites.

94

Actual locations
Estimated locations β = 0

12

x (cm)

10
8
6
4
2
0
0

1

2

1

2

3

4

5

6

7

3

4

5

6

7

5

6

7

y (cm)

3

2

1

0
0

Time (s)
(a)

Tracking error (cm)

0.4
0.3
0.2
0.1
0
0

1

2

3
4
Time (s)

(b)
Figure 4.16: Experimental results on the estimation of the source location: (a) Localization;
(b) tracking error.

95

Chapter 5
Underwater Tracking and
Size-Estimation of a Moving Object
In this chapter we investigate the problem of tracking a moving but non-vibrating cylindrical
object and estimating its size and shape using an artiﬁcial lateral line system. Based on a
nonlinear analytical model for the moving object-induced ﬂow ﬁeld, a two-stage extended
Kalman ﬁlter is proposed to estimate the location, velocity, size, and shape of the object.
For comparison purposes, unscented Kalman ﬁltering and particle ﬁltering are also explored.
The ﬂow model is validated with digital particle image velocimetry (DPIV) measurement.
Simulation results on tracking a cylinder with ellipsoidal cross-section are presented to illustrate the approach. The approach has been further veriﬁed experimentally with an artiﬁcial
lateral line prototype consisting of six ionic polymer metal composite (IPMC) ﬂow sensors.
The remainder of this chapter is organized as follows. We ﬁrst present the model for
the ﬂow ﬁeld generated by a moving object in a two-dimensional (2D) potential ﬂow, and
formulate the estimation problem in Section 5.1. The ﬁltering schemes are described in
96

Section 5.2, followed by the simulation example in Section 5.3 . The experimental setup and
results are presented in Section 5.4. Finally concluding remarks are provided in Section 5.5.

5.1

Flow Model and Problem Formulation

In this dissertation we assume a two-dimensional (2D) potential ﬂow. Consider a cylindrical
object moving through an otherwise still ﬂuid. We ﬁrst present the special case of circular
cross-section, and then generalize the model to an arbitrary cross-section proﬁle using the
conformal mapping theory.
Consider the 2D plane z = x+iy. Without the loss of generality, assume that the cylinder
is moving along x−direction with its cross-section lying in the x − y plane. The complex
potential wc (z), where z is outside the region occupied by the cylinder, is given by [77]

w c (z) = vx

R2
,
z − z1

(5.1)

where the superscript c indicates the case of a circular proﬁle, vx denotes the moving speed,
R is a radius of the circular cross-section, and z1 = xs + iys denotes the center of the moving
cylinder. The corresponding complex ﬂow velocity W c is then given by

W c (z) =

R2
dw c (z)
= −vx
.
dz
(z − z1 )2
97

(5.2)

c
c
If we write W c (z) = vx + i vy , then

vx R2 ((x − xs )2 − (y − ys )2 )
,
((x − xs )2 + (y − ys )2 )2
2vx R2 (y − ys )(x − xs )
c =
.
vy
((x − xs )2 + (y − ys )2 )2

c
vx = −

(5.3)
(5.4)

For a cylinder with a general shape (Fig. 5.1), its cross-section proﬁle can be obtained by
mapping the circular proﬁle with the Laurent series expansion [11]:

ζ(z) = (z − z1 ) +

λ1
λ2
+ ...,
+
(z − z1 ) (z − z1 )2

(5.5)

where λ1 , λ2 , · · · , are the shape parameters. When the shape parameters are real, the
obtained shape will be symmetric with respect to the x−direction. It can be shown that the
complex ﬂow velocity W g (z) around the moving object with a general proﬁle is

W g (z) =

dw c (z) dz
dz dζ

= vx −

R2
(z − z1 )2

W c (z)

1−

−1
λ1
− ...
.
(z − z1 )2

(5.6)

dz
dζ

In (5.6), z ∈ D, the domain exterior to the object with the general proﬁle. Since the impact
of higher-order terms in (5.5) decays quickly with the distance from the object [11], in this
chapter, we consider the case of ellipsoidal proﬁle only, namely, λ1 = 0, λi = 0, ∀ i > 1.
Furthermore, we assume that λ1 is real. Using (5.6), we can derive the complex ﬂow velocity

98

2

1

1
R

y (cm)

y (cm)

2

0
−1
−2
−2

0

−1

−1

0

1

−2
−2

2

x (cm)

−1

(a)

1

2

1

2

(b)
2

1

1

y (cm)

2

y (cm)

0

x (cm)

0

−1

−2
−2

0

−1

−1

0

1

−2
−2

2

x (cm)

−1

0

x (cm)

(c)

(d)

Figure 5.1: Illustration of the Laurent series expansion : (a) Circle ζ(z) = z; (b) the shape
ﬁngerprint ζ(z) = z + 1/(2z); (c) the shape ﬁngerprint ζ(z) = z + 1/(2z) + 1/(6z 2 ); (d) the
shape ﬁngerprint ζ(z) = z + 1/(2z) + 1/(6z 2 ) + 1/(12z 3 )(adapted from [11]).

W e (z) around an ellipsoidal cylinder:

W e (z) = vx −

R2
(z − z1 )2
99

1−

−1
λ1
.
(z − z1 )2

(5.7)

e
e
If we write W e (z) = vx + i vy , then

vx R2 (x − xs )2 − (y − ys )2 − λ1

e
vx = −
e
vy

=

(x − xs )2 − (y − ys )2 − λ1

2

+ 4(x − xs )2 (y − ys )2

2vx R2 (y − ys )(x − xs )

(x − xs )2 − (y − ys )2 − λ1

2

+ 4(x − xs )2 (y − ys )2

.

,

(5.8)
(5.9)

We now formulate the problem of localizing a moving object with ellipsoidal proﬁle and
estimating its size/shape in the x − y plane using an array of ﬂow sensors. Consider Fig. 5.2,
where an artiﬁcial lateral line comprising N sensors is located in parallel to the x-axis, with
the sensor locations denoted as (xi , yi ), 0 ≤ i ≤ N − 1. We denote the center of the moving
cylinder at time t by (xs (t), ys (t)). The object is assumed to move with a constant speed
vx parallel to the lateral line. We further assume that the presence of sensors has negligible
eﬀect on the ﬂow ﬁeld generated by the moving object. Each sensor is assumed to provide

(x (t ), y (t ))
1

y

r (t )
i

1

1

Moving object source

r (t )
i

2

(xi , yi )
Sensor i
x

Figure 5.2: Illustration of the problem setup for the tracking and estimation of a moving
object.
100

a noisy measurement of the local ﬂow velocity vx along the x−direction. In particular, the
measurement mi (t) by sensor i at time t is given by

mi (t) = f e (xi , yi , xs (t), ys (t), vx , R, λ1 ) + di (t),

(5.10)

where the function f e follows from (5.8), and di is the measurement noise. We introduce a
compact notation for the measurements from all sensors,

M (t) = [m1 (t) · · · mN (t)]T = F e (xs (t), ys (t), vx , R, λ1 ) + d(t),
where T denotes transpose,

△

F e (xs , ys , vx , R, λ1 ) =

[f e (x1 , y1 , xs , ys , vx , R, λ1 ) · · ·
f e (xN , yN , xs , ys , vx , R, λ1 )]T ,

and d(t) = [d1 (t) · · · dN (t)]T .
The estimation problem is then formulated as: given the measurements from the artiﬁcial
lateral line, M (·), determine the object location (xs (t), ys (t)), speed vx , size parameter R,
and shape parameter λ1 .

5.2

Nonlinear Filtering

Several ﬁltering schemes will be explored in this chapter, including extended Kalman ﬁltering,
unscented Kalman ﬁltering, and particle ﬁltering. We ﬁrst provide a brief review of each
method.
101

5.2.1

Extended Kalman Filter

In this chapter we propose to solve the tracking and estimation problem with extended
Kalman ﬁltering. We ﬁrst provide a brief review of extended Kalman ﬁltering, when the
measurement equation involves nonlinearities [76]. The discrete-time setting is considered
in this work, where the time index is denoted by k. Suppose that the system dynamics is
given by

Xk = Ak−1 Xk−1 + wk−1 ,
Mk = f (Xk ) + dk .

(5.11)

where Xk ∈ Rn denotes the state vector, Ak−1 has dimensions n × n, Mk ∈ Rp denotes the
measurement, f (·) is a nonlinear function, and wk−1 ∈ Rn and dk ∈ Rp denote the process
noise and the measurement noise, respectively. It is assumed that wk−1 and dk are white,
zero-mean, Gaussian noises, with covariance matrices denoted as Qk and Ck , respectively.
It is also assumed that wk−1 and dk are uncorrelated with each other.
Since f (Xk ) is a nonlinear function, in extended Kalman ﬁltering, it is linearized at
ˆ
X k|k−1 , the prediction for Xk at time k − 1. In particular, f (Xk ) is approximated by
ˆ
ˆ
ˆ
f (Xk ) ≈ f (X k|k−1 ) + Hk (Xk − X k|k−1 ),

(5.12)

ˆ
ˆ
where Hk is the Jacobin matrix of f (·) evaluated at X k|k−1 :
∂f (Xk )
ˆ
Hk =
.
∂Xk X =X
ˆ
k
k|k−1
102

(5.13)

ˆ
At k = 0, the initial state estimate X0|0 = X0 is a random vector with known mean
µ0 = E[X0 ], and the initial covariance matrix is given by P0|0 = E[(X0 − µ0 )(X0 − µ0 )T ].
The main steps in an extended Kalman ﬁlter are outlined next. The state prediction is given
by
ˆ
ˆ
Xk|k−1 = Ak−1 Xk−1|k−1 .

(5.14)

For k ≥ 1, the prediction of state covariance matrix follows
Pk|k−1 = Ak−1 Pk−1|k−1 AT + Qk−1 ,
k−1

(5.15)

and the optimal gain matrix Kk is given by

ˆ
ˆT
ˆT
Kk = Pk|k−1 Hk Ck + Hk Pk|k−1 Hk

−1

.

(5.16)

Finally, the estimates of the state and its covariance matrix are updated via

ˆ
Xk|k

ˆ
ˆ
= Xk|k−1 + Kk Mk − f Xk|k−1

Pk|k

ˆ
= (I − Kk Hk )Pk|k−1 ,

,

(5.17)
(5.18)

where I denotes the identity matrix.

5.2.2

Unscented Kalman Filter

Unlike the EKF, the Unscented Kalman Filter (UKF) does not approximate nonlinear function f (Xk ). Instead, it approximates the posterior p(Xk |Mk ) by a Gaussian density, which
is represented by a set of deterministically chosen sample points [76]. Consider the non103

linear ﬁltering problem deﬁned by (5.11). The assumption is that the posterior density at
ˆ
time k − 1 is Gaussian: p(Xk−1 |Mk−1 ) = N (Xk−1 ; Xk−1|k−1 , Pk−1|k−1 ). The ﬁrst step is
i
i
to represent this density by a set of Ns sample points Xk−1|k−1 and their weights Wk−1 ,

i = 0, . . . , Ns − 1. The prediction step is then performed as follows:
ˆ
Xk|k−1 =

Ns −1
i=0

i
i
Wk−1 · Xk−1|k−1 ,

Pk|k−1 = Qk−1 +

Ns −1
i=0

(5.19)

i
i
ˆi
Wk−1 Xk−1|k−1 − Xk|k−1

i
ˆi
Xk−1|k−1 − Xk|k−1

T

.
(5.20)

ˆ
The predicted density p(Xk |Mk−1 ) ≈ N (Xk ; Xk|k−1 , Pk|k−1 ) is represented by a set of Ns
sample points, given as
i
i
Xk|k−1 = Ak−1 Xk−1|k−1 ,

(5.21)

the predicted measurement step is given by

ˆ
Mk|k−1 =

Ns −1
i=0

i
i
Wk−1 fk Xk|k−1 ,

(5.22)

and the update step is given as

ˆ
Xk|k

ˆ
ˆ
= Xk|k−1 + Kk Mk − Mk|k−1 ,

(5.23)

Pk|k

T
= Pk|k−1 − Kk Sk Kk ,

(5.24)

104

where


Ns −1

Kk = 

i=0

Sk = C k +

i
i
ˆ
Wk−1 Xk|k−1 − Xk|k−1

Ns −1
i=0

i
ˆ
fk Xk|k−1 − Mk|k−1

i
i
ˆ
Wk−1 fk Xk|k−1 − Mk|k−1

T



 S −1 ,

(5.25)

k

i
ˆ
fk Xk|k−1 − Mk|k−1

T

.
(5.26)

ˆ
The n dimensional random variable X(k) with mean X (k |k ) and covariance P (k |k ) is
approximated by 2n weighted samples or sigma points selected by the algorithm [76]

α
W0 = (n+α) ,

ˆ
X0 (k |k ) = X (k |k )
ˆ
Xi (k |k ) = X (k |k ) +
ˆ
Xi+n (k |k ) = X (k |k ) −

(n + α)P (k |k )
(n + α)P (k |k )

where α is scaling parameter, such that α + n = 0,

i

α
Wi = 2(n+α) ,
α
Wi+n = 2(n+α) ,
i

(n + α)P (k |k )

i

is the ith column of

the lower Cholesky factorization of (n + α)P (k |k ), and Wi is the weight that is associated
with the ith point.

5.2.3

Particle Filter

Particle ﬁlter (PF) is a suboptimal ﬁlter. It performs sequential Monte Carlo (SMC) estimation based on particle representation of probability densities. In this chapter, we have
chosen the modiﬁed particle ﬁlter, called regularized particle ﬁlter (RBF) [62]. Regularization consists of resampling the particles according to a continuous approximation of the
target distribution so that all the particles obtained have diﬀerent locations. The smoothing
is performed by convoluting the discrete PF approximation with a kernel whose properties
105

are provided by the density estimation theory [85]. The resulting continuous approximation
takes the following form:
Ns

p (Xk |M1:k ) ≈
ˆ

i=1

i
i
wk Kβ (Xk − Xk ),

(5.27)

where
1
Kβ (X) = n K
β

X
β

(5.28)

is the rescaled kernel function density K(·), β > 0 is the kernel bandwidth (a scalar parami
eter), Ns is the total number of particles, and wk , i = 0, . . . , Ns are the normalized weights.

The particle system evolves according to an importance-sampling rule; more precisely, the
particles are simulated sequentially according to a proposal distribution given as

i
i
Xk ∼ q Xk X1:k−1 , M1:k

f or i = 1, . . . , Ns ,

(5.29)

where q(·) is a proposal density. They are then assigned importance weights to correct for
the discrepancy between the proposal and the target distributions given by
Ns

wk = w i /
˜

wj ,
˜

(5.30)

j=1

where
wi =
˜

i
i
i
i
wk−1 p Mk X1:k , M1:k−1 p Xk X1:k−1
.
q (Xk |X1:k−1 , M1:k )

(5.31)

According to (5.31), the most likely particles yield high importance weights. A selection
step is ﬁnally introduced to prevent degeneracy. This selection is performed by resampling
the set of particles according to the obtained approximation of the posterior distribution.
106

However, systematic resampling is known to result in a loss of sample diversity. A measure
of degeneracy, called eﬀective sample size and denoted Nef f , has been introduced in [51] to
decide whether resampling is useful or not:

Nef f =

1
Ns
i 2
i=1 wk

(5.32)

The selection procedure is carried out whenever Nef f is below a given threshold.
For our tracking and estimation problem, the dynamics of interest can be described as
follows:

xs [k] = xs [k − 1] + vx [k − 1]∆ + w1 [k − 1],

(5.33)

ys [k] = ys [k − 1] + vy [k − 1]∆ + w2 [k − 1],

(5.34)

vx [k] = vx [k − 1] + w3 [k − 1],

(5.35)

vy [k] = vy [k − 1] + w4 [k − 1],

(5.36)

R[k] = R[k − 1] + w5 [k − 1],

(5.37)

λ1 [k] = λ1 [k − 1] + w6 [k − 1],

(5.38)

where ∆ is the sampling time, and w1 , · · · w6 denote the process noises, which are assumed
to be uncorrelated. The measurement equation is given by

M [k] = F e (xs [k], ys [k], vx [k], R[k], λ1 [k]) + d[k],

(5.39)

where d[k] denotes the vector of measurement noises with uncorrelated components.
In this work, we propose a two-stage nonlinear ﬁltering scheme, as illustrated in Fig. 5.3.
107

In the ﬁrst stage, ﬁve state components, (5.33) – (5.37), are estimated while assuming the
shape parameter λ1 = 0. Correspondingly, the measurement equation for the ﬁrst stage is
obtained by plugging λ1 = 0 in (5.39):

M [k] = F e (xs [k], ys [k], vx [k], R[k], 0) + d[k].

(5.40)

In the second stage, the state estimate involves λ1 only (equation (5.38)), while the measurement equation is obtained from (5.39) by plugging in the state estimates from the ﬁrst
stage:
ˆ
M [k] = F e (ˆs,k|k , ys,k|k , vx,k|k , Rk|k , λ1 ) + d[k].
x
ˆ
ˆ

(5.41)

State estimation at each stage is then carried out following the general nonlinear ﬁltering
procedure, as outlined earlier.

Figure 5.3: Illustration of the two-stage ﬁltering scheme.
There are two motivations for estimating the desired information in two cascaded stages.
108

First, as discussed in [11], it is likely that a biological lateral line extracts the relevant
information in a progressive manner. Second, the proposed scheme leads to reduced computational complexity because of the reduced dimensions for the state dynamics – for example,
one only needs to invert 5 × 5 matrices instead of 6 × 6 matrices.

5.3
5.3.1

Simulation Results
Recursive Information Matrix

Generally, a nonlinear ﬁltering problem does not have a close-form analytical solution, and
in all practical applications, nonlinear ﬁltering is performed by some sort of approximation.
Despite the absence of a closed-form solution, we would like to ﬁnd the best achievable error
performance for nonlinear ﬁltering; to do so, a theoretical Cramer-Rao lower bound has been
formulated for such problems. Consider the system which is described in (5.11), the recursive
information matrix Jk for such system has been derived in [76] as

˜ T −1 ˜
Jk = Q−1 + Hk Ck Hk − Q−1 Ak−1 Jk−1 + AT Q−1 Ak−1
k−1 k−1
k−1
k−1

−1

AT Q−1 ,
k−1 k−1

(5.42)

where
∂f (Xk )
˜
H=
.
∂Xk

(5.43)

is the Jacobian of f (·) evaluated at the true value of Xk .
Fig. 5.4 shows the square root of the CRB,

√

CRB, corresponding to the components of

the state vector: (i) xs -coordinate, (ii) ys -coordinate, (iii) size parameter R , and (iv) shape
parameter λ1 , for two-stage and one-stage schemes, respectively. These bounds are obtained
109

as:
CRB(Xk [i]) = J −1 [i, i],

i = 1, 2, 3, 4.

(5.44)

The close-form analytical solution for our problem has derived in Appendix B.1. It is clear

(a)

(b)

4

4
Two stage
One stage

2

2

y

s

σ (cm)

3

s

σx (cm)

3

Two stage
One stage

1
0
0

1

0.5

1

0
0

1.5

0.5

Time (s)
(c)

1.5

(d)

5

5
Two stage
One stage

Two stage
One stage

4

σλ (cm)

4
3

3

1

σR (cm)

1

Time (s)

2

1

1
0
0

2

0.5

1

0
0

1.5

0.5

1

1.5

Time (s)

Time (s)

√
Figure 5.4: CRB for: (a) Moving object xs coordinate; (b) moving object ys coordinate;
(c) size parameter R; (d) shape parameter λ1 .

from Fig. 5.4 that, a two-stage scheme yields smaller estimation errors than a one-stage
scheme. Consequently, we will employ two-stage schemes in our estimation.
110

5.3.2

Simulation Results on Nonlinear Filtering

Fig. 5.5 illustrates the simulation setup. The lateral line system is placed parallel to the
x−axis and centered at (5, 0) cm. It consists of 6 sensors, with the sensor-to-sensor separation
of 2 cm. The ellipsoidal cylinder, with a shape parameter (λ1 = 1.5 cm) and a size parameter
R = 2.5 cm, moves from left to right with a constant speed of 6 cm/s. The initial location of
the moving object is (0, 7) cm, and the time duration for the movement is 1.667 s, resulting
in a terminal location of (10, 7) cm. The initial state estimates are set to be xs,0|0 = 0.1 cm,
ˆ
ˆ
ˆ
ys,0|0 = 5 cm, vx,0|0 = 4 cm/s, vy,0|0 = 1 cm/s, R0|0 = 1 cm, and λ1,0|0 = 0.1 cm. The
ˆ
ˆ
ˆ
sampling time ∆ is chosen to be 0.01 s. Fig. 5.6 shows the simulated measurements of
ﬂow velocities (x−component) at all sensor sites, where a zero-mean, Gaussian noise with a
2
variance of (σm = 0.01 cm2 /s2 ) is added to each measurement.

9
8 (x (t ), y (t ))
s 1

7

(xs(tn ), ys(tn ))

(x (t ), y (t ))
s f

s f

Vx

6

y (cm)

s 1

5
4
3
2

Sensors

1
0

0

2

4

6

8

10

12

x (cm)

Figure 5.5: Simulation setup for tracking of a moving cylinder.
Fig. 5.7(a) shows the actual object locations and those estimated using the proposed
two-stage extended Kalman ﬁltering scheme for a typical run. Here three diﬀerent val111

Sensor#1
Sensor#2
Sensor#3
Sensor#4
Sensor#5
Sensor#6

Velocity (cm/s)

0.4
0.3
0.2
0.1
0
−0.1
0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

Time (s)

Figure 5.6: Simulated time trajectories of sensor measurements.
ues of variance σ 2 for the process noises (w1 , · · · , w6 ) are tested, to examine the impact
of the process noise variance on the estimation results. Fig. 5.7(b) shows the localization error for each actual x−coordinate (xs [k]) of the object, where the error is deﬁned as
(ˆs,k|k − xs [k])2 + (ˆs,k|k − ys [k])2 .
x
y
It can be observed that the estimated object location quickly converges to the neighborhood of the actual location after a short transient period (less than 0.15 s). In addition, while
the variance of the process noises does inﬂuence the transients, its impact on the estimation
performance at the steady state is not pronounced. The insensitivity of the estimation performance to the process noise variances is a positive news, since it is diﬃcult to precisely
characterize these noise variance values in practice. Table 5.1 summarizes the averages of
the estimates for other parameters (vx , R, and λ1 ), where each average is obtained by taking
112

the mean of the estimated values over the period [0.33, 1.67] s.

We have further studied

Actual Track
Estimated track (σ2 = 0.1)
Estimated track (σ2 = 0.3)
Estimated track (σ2 = 0.5)

8

y (cm)

6
4
2
0
0

2

4

6

8

10

x (cm)
(a)

Localization error (cm)

1
2

σ = 0.1

0.8

2

σ = 0.3
2

σ = 0.5

0.6
0.4
0.2
0
0

2

4

x (cm) 6

8

10

(b)
Figure 5.7: Simulation results on tracking the ellipsoidal cylinder using two-stage extended
Kalman ﬁltering: (a) The estimation of the moving object location for diﬀerent variance
values for the process noises; (b) localization error at each x-coordinate of the object. The
measurement noise variance is held at 0.01 cm2 /s2 .
2
the impact of the measurement noise variance σm . Fig. 5.8 shows the localization results for
2
diﬀerent values of σm , while the variance for each process noise, σ 2 , is held constant at 0.3

113

Table 5.1: Simulation results: Estimated velocity, size parameter, and shape parameter
under diﬀerent process noise variances. Two-stage extended Kalman ﬁltering is used.

Parameter
Velocity (cm/s)
Size R (cm)
Shape λ1 (cm)

Actual
6
2.5
1.5

Estimated (σ 2 = 0.1)
6.22
2.34
1.23

(σ 2 = 0.3)
6.35
2.41
1.42

(σ 2 = 0.5)
6.44
2.32
1.43

Table 5.2: Simulation results: Estimated velocity, size parameter, and shape parameter
under diﬀerent measurement noise variances. Estimation is done with two-stage extended
Kalman ﬁltering.

Parameter
Velocity( cm/s)
Size R (cm)
Shape λ1 (cm)

Actual
6
2.5
1.5

2
Estimated (σm = 0.01)
6.45
2.38
1.45

2
(σm = 0.03)
7.23
2.24
1.42

2
(σm = 0.05)
7.24
2.258
1.63

(in its respective unit), while Table 5.2 lists the averages of estimated values for vx , R, and
λ1 .
It can be seen that, overall, the measurement noise variance has an appreciable impact on
the estimation performance; in general, larger measurement noise variance leads to longer
transient and larger estimation error, which is consistent with one’s intuition. We have
further studied the impact of the input and measurement noise variances for a one-stage
extended Kalman ﬁltering scheme in Appendix B.2.
Finally, the performance of the three nonlinear tracking ﬁlters, designed for this problem,
are compared against each other. The considered ﬁlters include the extended Kalman ﬁltering
(EKF), unscented Kalman ﬁltering (UKF), and regularized particle ﬁltering (RPF). Under
two-stage typical run, the results for the actual object locations and those mean estimated
and their estimation errors are shown in Fig. 5.9.
114

From Fig. 5.9, we can see that the three nonlinear ﬁlters are unbiased and all of them
demonstrate almost similar tracking accuracy. Table 5.3 summarizes the averages of the
estimates for other parameters (vx , R, and λ1 ) for the three nonlinear ﬁlters.
Table 5.3: Simulation results: Estimated velocity, size parameter, and shape parameter
under three nonlinear ﬁltering algorithms.
Parameter
Velocity (cm/s)
Size R (cm)
Shape λ1 (cm)

Actual
6
2.5
1.5

Estimated (EKF)
6.45
2.83
1.45

Estimated (UKF)
6.33
2.24
1.42

Estimated (RPF)
6.23
2.25
1.55

In terms of the computational load of the three algorithms, the following results are
obtained (based on MATLAB implementation). The UKF requires 5 times more CPU time
than the EKF. The RPF (with 1000 particles) requires 110 times more time than the EKF.
Therefore, the EKF is the most eﬃcient scheme among those compared.

5.4

Experimental Results

Fig. 5.10 shows the schematic of the experimental setup, while Figs. 5.11(a) and (b) show
the picture of the actual experimental system. The artiﬁcial lateral line consisted of six
IPMC sensors. The sensor-to-sensor separation was 2 cm, resulting in a total span of 10 cm.
Each sensor was mounted in such a way that it would bend in the direction parallel to the
sensor array, which is denoted as x−direction.

The experiments were conducted in a water

tank measuring 6 ft long, 2 ft wide, and 2 ft deep. A digital particle image velocimetry
(DPIV) system (LaVision, Ypsilanti, MI), not shown in Fig. ??, which was explained earlier
in Chapter 3, was used to measure the ﬂow for model validation. A custom-built, motordriven conveyor was used to move an object in the water. The system can be driven by a DC
115

motor or an AC motor. In our experiments, two objects were used, the ﬁrst one is a circular
cylinder with outer radius of 3.02 cm and an ellipsoidal cylinder with the minor and major
diameters of 2.5 cm and 5 cm, respectively. Both objects had a height of 10 cm. The speed
of the object can be adjusted by varying the voltage to the motor or changing the setup of
the gearbox that is associated with the AC motor. The conveyor track was separated from
the lateral line system by 7 cm. The signals from the IPMC sensors were conditioned and
ampliﬁed before the acquisition and processing by a PC equipped with dSPACE (DS 1104,
dSPACE).
We ﬁrst conducted DPIV measurements (using DC motor driving the moving object) to
validate the ﬂow model (5.3). Fig. 5.12 shows the comparison of the measured ﬂow velocity
(x−component) with the prediction by (5.3) at a particular sensor site, when the circular
cylinder was moved to pass by the sensor site at a speed of 6.67 cm/s. The reasonably good
match between the model and the measurement in Fig. 5.12 supports the use of the ﬂow
model in this work.
Before running the tracking experiments, we performed the calibration of all IPMC sensors using the mini-shaker. The details of the calibration process are explained earlier in
Chapter 3.
Fig. 5.13 shows the measured ﬂow velocities by all six IPMC sensors when the cylindrical and ellipsoidal object traveled by the lateral line system at 6.67 cm/s and 40 cm/s,
respectively, which demonstrates a similar shape pattern as in Fig. 5.6.
Fig. 5.14(a) shows the actual trajectory of the moving cylinder and its estimate using the
ﬁrst-stage Kalman ﬁltering shown in Fig. 5.3, and Fig. 5.14(b) shows the tracking error. Since
it was known that the cross-section of this cylinder was circular, the shape parameter was
116

not identiﬁed in these experiments. The process noise variances were chosen to be σ 2 = 0.3
(in its respective unit of each state component) and the measurement noise variances were
2
set to be σm = 0.01 (cm/s)2 . Table 5.4 lists the estimated moving speed and size parameter.

The tracking error was less than 0.5 cm and the errors in estimating vx and R were less
than 25%, which demonstrates the eﬀectiveness of the proposed modeling and estimation
approach.
Table 5.4: Cylindrical object experimental results: Estimated velocity and size parameter.
Parameter
Velocity (cm/s)
Size R (cm)

Actual
6.667
3.02

Estimated
8.3
3.51

Error
1.6
0.49

Fig. 5.15 (a) shows the actual trajectory of the moving ellipsoidal object and its estimate
using the two-stage Kalman ﬁltering shown in Fig. 5.3, and Fig. 5.14(b) shows the tracking
error. The actual size and the shape parameter of the ellipsoidal object were identiﬁed to be
R = 1.5 cm and λ1 = 1.25 cm, respectively. The process noise variances and measurement
noise were chosen as explained earlier in the cylindrical experiment setup.
Table 5.5 lists the estimated moving speed, size, and shape parameter. The tracking
error was less than 0.45 cm and the errors in estimating vx and R were less than 33%, which
demonstrates the eﬀectiveness of the proposed modeling and estimation approach.
Table 5.5: Ellipsoidal object experimental results: Estimated velocity, size parameter, and
shape parameter.
Parameter
Velocity (cm/s)
Size R (cm)
Shape λ1 (cm)

Actual
40
1.5
1.25

117

Estimated
41
1.9
1.6

Error
1
0.4
0.35

5.5

Chapter Summary

In this chapter we proposed new nonlinear ﬁltering algorithms for object tracking and estimation using artiﬁcial lateral lines. Based on a 2D potential ﬂow model, an estimation
problem was proposed and solved with a two-stage nonlinear ﬁltering method. Extended
Kalman ﬁltering was shown to have similar tracking accuracy performance as unscented
Kalman ﬁltering and regularized particle ﬁltering while having the least computational complexity. Simulation results were presented to illustrate and demonstrate the eﬀective of the
proposed approach, where the impacts of process noises and measurement noises were explored. Experimental results were also reported, where an artiﬁcial lateral line prototype
consisting of six IPMC sensors was used to track and estimate the size and the shape of both
a cylinder and an ellipse moving object.

118

Actual track
2

Estimated track (σm = 0.01)
2

Estimated track (σm = 0.03)
2

Estimated track (σm = 0.05)

y (cm)

8
6
4
2
0
0

2

4

6

8

10

x (cm)
(a)
2

σm = 0.01

Localization error (cm)

2

σm = 0.03

1

σ2 = 0.05
m

0.8
0.6
0.4
0.2
0
0

2

4

x (cm)

6

8

10

(b)
Figure 5.8: Simulation results on tracking the ellipsoidal cylinder using the two-stage extended Kalman ﬁltering: (a) The estimation of the moving object location for diﬀerent
variance values for the measurement noises; (b) localization error at each x-coordinate of
the object. The process noise variance is held at 0.3 (in its respective unit for each state
component).

119

Actual track
EKF
UKF
RPF

8

y (cm)

6
4

Localization error (cm)

2
0
0
1

2

4

2

4

6

8

10

6

8

10

0.8
0.6
0.4
0.2
0
0

x (cm)

Figure 5.9: Simulation results on tracking the ellipsoidal cylinder for three nonlinear ﬁlters:
The process noise variance is held at 0.3 (in its respective unit for each state component)
and the measurement noise variance is held at 0.01 cm2 /s2 .

120

Power
supply
dSPACE
Cylinderical
object
DC Motor driving
the moving object

Electronic
conditioning

Supporting
frame
rod
IPMC Sensors
IPMC bending
direction
Moving direction
Tank

Figure 5.10: Experimental setup: Schematic.

121

Custom built
conveyor

DC Motor

Cylinderical
object

Supporting
frame

IPMC sensors

(a)

AC Motor

Ellipsoidal
object
Custom built
conveyor
Supporting
frame

IPMC
sensors

(b)
Figure 5.11: Experimental setup: (a) Picture of the experimental setup with a circular
cylinder; (b) picture of the experimental setup with an ellipsoidal cylinder.

122

1.5
Measurment
Model

Velocity (cm/s)

1
0.5
0
−0.5
−1
0

Sensor

0.5

1

1.5

2

2.5

3

3.5

4

Time (s)
Figure 5.12: Comparison of DPIV-measured ﬂow velocity with the model prediction.

123

Sensor #1
Sensor #2
Sensor #3
Sensor #4
Sensor #5
Sensor #6

Velocity (cm/s)

2
1.5
1
0.5
0
0

0.5

Time (s)

1

1.5

(a)

Sensor #1
Sensor #2
Sensor #3
Sensor #4
Sensor #5
Sensor #6

Velocity (cm/s)

3
2
1
0
−1
0

0.05

0.1

0.15

0.2

0.25

0.3

Time (s)
(b)
Figure 5.13: Measured ﬂow velocities by six IPMC sensors: (a) DC motor driving the
cylindrical object, with digital low pass ﬁlter (fc = 42) Hz; (b) AC motor driving ellipsoidal
object, with digital low pass ﬁlter (fc = 4) Hz.
124

Actual track
Estimated track

8

y (cm)

6
4
2
0
0

2

4

x (cm)

6

8

10

6

8

10

Localization error (cm)

(a)

1
0.8
0.6
0.4
0.2
0
0

2

4

x (cm)
(b)

Figure 5.14: Experimental results on tracking a cylinder with circular proﬁle with two-stage
extended Kalman ﬁltering: (a) The estimation of the object location; (b) localization error
at each x−coordinate of the object.

125

Actual track
Estimated track

y (cm)

8
6
4
2
0
0

2

4

6

8

10

6

x (cm)

8

10

Localization error (cm)

(a)

1
0.8
0.6
0.4
0.2
0
0

2

4

x (cm)
(b)
Figure 5.15: Ellipsoidal object experimental results on tracking an ellipse: (a) The estimation
of the object location; (b) localization error at each x−coordinate of the object.

126

Chapter 6

Exploration of IPMC Sensors in
Vortex Street Sensing

Most ﬁsh frequently exploit unsteady forcing of ambient/wake ﬂow to extract energy for
eﬃcient locomotion and for obstacle avoidance. In this chapter, we ﬁrst describe an analytical
model for a pair of stable vortices outside an obstacle and then for a K´rm´n vortex street.
a a
Second, with extensive experiments, we explore the problem of identifying hydrodynamic
parameters associated with on obstacle in a ﬂow using IPMC sensors. In particular, we show
that it is feasible to identify ﬂow velocity and obstacle size based on the measurement of
vortex shedding frequency.
The remainder of the chapter is organized as follows. In Section 6.1 the model for a
stationary vortex and that for a K´rm´n vortex street are described. The localization of a
a a
static obstacle is discussed in Section 6.2. The experimental setup and results are described
in Section 6.3. Finally, chapter summary is provided in Section 6.4.
127

6.1

Stationary Vortex and K´rm´n Vortex Street
a a

The circle theory reported in [89] proves that any pair of vortices behind a cylinder is stable.
Therefore, for every point z outside the circle of radius R in the complex plane, there is
another opposite point (its complex conjugate), z = R2 /z, which lies inside the circle. Thus,
¯
we have two singularities at these points: f (z), where a point z lies outside the circle, and
¯z
¯
f (¯) = f (R2 /z), where a point z lies inside the circle. The complex potential w(z) of these
¯
singularities is given by
¯
w(z) = f (z) + f (R2 /z).

(6.1)

Generally, the circle theory says that for each vortex outside the circle has a vortex of equal
strength but opposite direction inside the circle, which provides an analytical tool for modeling vertex shedding behind a circular cylindrical obstacle. First, we consider two vortices
located in the complex plane. The ﬁrst vortex at position p1 has a complex coordinate z1
and strength k, and the other vortex at position p2 has a complex coordinate z1 and strength
¯
−k, as shown in Fig. 6.1.

The complex potential w in presence of the cylinder |z| = R is

given as

w = ik ln (z − z1 ) + wz ,
wz = −ik ln (z − z1 ) − V
¯

z + R2 /z

− ik ln

R2 − z
¯
1
z
.
R2 − z
1
z

(6.2)

Here V is the nominal cross ﬂow velocity. The vortex at p1 will be at rest if dwz = 0, which
dz
means that its velocity is zero. With this condition, it can be demonstrated [89] that the
128

y

p

1

p

1

o

x

p2
p

2

Figure 6.1: Illustration of the analytical model for cylinder-shed vorticity.

vortices behind a cylinder in ﬂow are stable with the vortex strength of k given by

k = V ( r 2 − R2 )( r 2 + R2 )/ r 5 ,

(6.3)

where r = (x, y)T , and T denotes the transpose. Accordingly, we can ﬁnd the exact position
of the stable vortices behind the cylinder when we use the transformation z = r eiθ , where
0 ≤ θ ≤ 90o , implying
r 2=

R2
.
1 − sin(θ)

(6.4)

Fig. 6.2 shows the stable vortices behind a cylinder (R = 1 cm) in a ﬂow with nominal
velocity of 2 cm/s.
Second, we consider the case of a shed K´rm´n vortex street which consists of two parallel
a a
129

y (cm)

5

0

−5
−5

0
x (cm)

5

Figure 6.2: Visualization of the ﬂow around a cylinder.
semi inﬁnite columns of vortices of the same spacing, say a. We denote the vertical distance
between the columns as b. The vortex strength is k for all vortices with y > 0 and it is −k
for all vortices with y < 0 as shown in Fig. 6.3.

y
a

V Flow

R

b/2

x
a

Figure 6.3: Illustration of the analytical model for cylinder-shed K´rm´n vortex street.
a a
130

The complex potential for the vortex street at instant t = 0 is given as [89],

c
wo (z) = V (z + R2 /z) + ik ln sin

π
a

z−

ib
2

c
which implies that complex ﬂow velocity Wo =

c
Wo (z) =

V (1 − (R/z)2 ) +

iπk
a

cot

π
a

− ik ln sin

π
a

z−

a ib
+
2
2

,

(6.5)

c
dwo (z)
at time t = 0 is given by
dz

z−

ib
2

− cot

π
a

z−

a ib
+
2
2

, (6.6)

where z is outside the region occupied by the cylinder. The complex ﬂow velocity at any
instant t is given by

c
Wo (z) = V (1 − (R/z)2 ) +

iπk
a

cot

π
a

z − uv t −

ib
2

− cot

π
a

z − uv t −

a ib
+
2
2

,(6.7)

where uv , the velocity of the vortex street, is given by

uv = V −

πk
tanh
a

b
π .
a

(6.8)

Finally, the relationship between a and b is approximately [89]
√
1
b
= cosh−1 2.
a
π

(6.9)

The initial conﬁguration will repeat itself after it has moved a horizontal distance equal to
a.

Now, we need to relate this solution to an ellipsoidal obstacle shape. To accomplish this,
we use a conformal mapping function called the Joukowsky transformation. The Joukowsky
131

transformation transforms circles in the ζ plane into shapes that resemble ellipses in the z
plane [77]. In contrast to the usual way, we will not study the ﬂow velocity in the ζ-plane,
but in the original z-plane. Consider, ζ(z) := z + λ2 /z with a real shape parameter λ, which
implies that the obtained obstacle shape will be symmetric with respect to the x-direction.
e
It can be shown that the complex ﬂow velocity Wo around an ellipsoidal obstacle is given

by

e
Wo (z) =

dζ −1
dz

dwc (z)
dz
c
Wo

=

V (1 − (R/z)2 ) +
1 − (λ/z)2

−1

iπk
a

cot

π
a

z − uv t −

ib
2

− cot

π
a

z − uv t −

.

a ib
+
2
2
(6.10)

In (6.10), z ∈ D, the domain outside an ellipsoidal obstacle. Fig. 6.4 shows the visualization
of a K´rm´n vortex street behind an obstacle, where Table 6.1 summarizes the values of the
a a
setup parameters.

Table 6.1: Simulation results: Setup parameters.
Parameter

Value

Nominal ﬂow velocity (cm/s)
Distance between vortices a (cm)
Vortex strength
Radius R (cm)
Shape parameter λ (cm)

2
1
1
1
0.5

132

2

y (cm)

4

2

y (cm)

4

0

0
−2

−2
−4
−5

0

−4
−5

5

x (cm)

(a)

0

x (cm)

5

(b)

Figure 6.4: Visualization of a K´rm´n vortex street ﬂow ﬁled around: (a) Circular cylinder;
a a
(b) ellipsoidal cylinder.

6.2

Localization of an Obstacle: Problem Formulation

We now deﬁned the problem in the two-dimensional x − y plane. As illustrated in Fig. 6.5,
suppose that a cylindrical static obstacle is located at (xo , yo ), where (zo = xo + iyo ) is
the center of that obstacle and its size parameter (radius of circle) is Ro , in a ﬂow with
constant nominal ﬂow velocity V that has an angle α with respect to the x-axis. It will
generate a vortex street with the vortex sheding frequency fsh . The Strouhal number (St)
is a dimensionless index, given by
f (2R)
St = sh c ,
v

(6.11)

where v c is the actual ﬂow velocity in the region behind the cylindrical obstacle. The Strouhal
number depends on the Reynolds number; however, it remains at an almost constant value
of 0.22 for Reynolds number > 2000 [4].

Consider Fig. 6.5, where an artiﬁcial lateral line

comprising N sensors is located in parallel to the x-axis, with the sensor locations denoted
133

Obstacle center
(x , y )

y

o

o

Sensors
(x , y )

Ro

i

V

i

Vortex

α

x

Figure 6.5: Illustration of the problem setup. The obstacle is located at (xo , yo ).
as (xi , yi ), 0 ≤ i ≤ N − 1. We assume that the presence of sensors has negligible eﬀect
on the ﬂow distribution. For a circular cylinder, we can derive the complex ﬂow velocity
c
c
c
Wo = vx + ivy around the cylinder:
2
2
2Ro (xi − xo )2
Ro
−
(xi − xo )2 + (yi − yo )2 ((xi − xo )2 + (yi − yo )2 )2
iπk
π
ib
+Re
cot
(xi − xo ) + i(yi − yo ) − uv t −
a
a
2
π
a ib
− cot
(xi − xo ) + i(yi − yo ) − uv t − +
,
(6.12)
a
2
2

c
vx (xi , yi , t) = V cos(α) 1 +

c
vy (xi , yi , t) = −V sin(α)

2
2Ro (xi − xo )(yi − yo )

2

(xi − xo )2 + (yi − yo )2
π
ib
iπk
cot
(xi − xo ) + i(yi − yo ) − uv t −
+Img
a
a
2
a ib
π
(xi − xo ) + i(yi − yo ) − uv t − +
.
− cot
a
2
2
134

(6.13)

e
e
e
Using Joukowski transformation, we can obtain the ﬂow velocity Wo = vx + ivy around an

ellipsoidal cylinder:

e
vx (xi , yi , t) = V cos(α) 1 +

−

2
(Ro − λ2 )((xi − xo ) + λ)
2λ ((xi − xo ) + λ)2 + (yi − yo )2

(R2 − λ2 )((xi − xo ) − λ)
2λ ((xi − xo ) − λ)2 + (yi − yo )2

2πk(xi − xo )(yi − yo )λ2

+

(xi − xo )2 − (yi − yo )2 − λ2

a
Re

cot

− cot

π
a

2

+ 4(xi − xo )2 (yi − yo )2

ib
2
a ib
(xi − xo ) + i(yi − yo ) − uv t − +
2
2

e
vy (xi , yi , t) = −V sin(α)

π
a

(xi − xo ) + i(yi − yo ) − uv t −

,

2
2(Ro − λ2 )(xi − xo )(yi − yo )

λ4 − 2λ2 (xi − xo )2 − (yi − yo )2 + (xi − xo )2 − (yi − yo )2

+

πk (xi − xo )2 + (yi − yo )2
a

π
a

2

+ 4(xi − xo )2 (yi − yo )2

ib
2
a ib
(xi − xo ) + i(yi − yo ) − uv t − +
2
2

cot

π
a

2

− λ2 (xi − xo )2 − (yi − yo )2

(xi − xo )2 − (yi − yo )2 − λ2

Img
− cot

2

(6.14)

(xi − xo ) + i(yi − yo ) − uv t −

.

(6.15)

Fig 6.6 shows a typical ﬂow velocity behind a circular cylinder (Ro =2 cm) and an ellipsoidal
obstacle shape (λ = 1.5 cm), respectively. The obstacle’s center is located at (xo = −2, yo =
5) cm. The obstacles are exposed to a cross ﬂow with nominal ﬂow velocity of 15 cm/s, and
the vortex strength is set to be k = 0.08.

Each sensor is assumed to provide a noisy measurement of the local ﬂow velocity vx along
135

20
ve (cm/s)

vx (cm/s)

20
10

x

c

0

0

−10
10
10

y (cm) 5

10

y (cm)5

5

0 0

−20
10

x (cm)

5
x (cm)

0 0

(a)

20

5

10

ve (cm/s)
y

vc (cm/s)
y

10

0
−5
−10
10

10

y (cm)

5

−10
−20
10
10
5

5
0 0

0

y (cm)

x (cm)

5
0 0

x (cm)

(b)
Figure 6.6: Simulation results: 3D illustration of ﬂow velocity behind the obstacles: (a)
Along the x−direction; (b) along the y−direction. In each subﬁgure, the left plot shows the
velocity for the circular cylinder while the right plot shows the velocity for the ellipsoidal
cylinder.

the x−direction. In particular, the measurement mi (t) by sensor i at time t is given by

mi (t) = f e (xi , yi , xo , yo , V, α, Ro , λ, a, b, uv , k) + di (t),
136

(6.16)

where the function f e follows from (6.14), and di is the measurement noise. We introduce a
compact notation for the measurements from all sensors,

M (t) = [m1 (t) · · · mN (t)]T = F e (xo , yo , V, α, Ro , λ, a, b, uv , k) + d(t),
where T denotes transpose,

△

F e (xo , yo , V, α, Ro , λ, a, b, uv , k) =

[f e (x1 , y1 , xo , yo , V, α, Ro , λ, a, b, uv , k) · · ·
f e (xN , yN , xo , yo , V, α, Ro , λ, a, b, uv , k)]T ,

and d(t) = [d1 (t) · · · dN (t)]T .
The estimation problem is then formulated as: given the measurements from the artiﬁcial
lateral line, M (·), estimate the obstacle location (xo , yo ), size Ro and shape parameter λ,
ﬂow velocity V , and angle α. In the next section, we provide some experimental evidence
on how part of the estimation problem can be solved.

6.3
6.3.1

Experimental Setup and Results
Swim Tunnel Flow Calibration

We have conducted experiments in a swim tunnel (Loligo Systems) measuring (length = 146
× width = 68 × height = 35) cm3 , with a working section of (65 × 20 × 20) cm3 , as shown
in Fig. 6.7.
The nominal ﬂow velocity in the working section of the swim tunnel tank is controlling
through a knob. In order to calibrate the swim tunnel, we measured the nominal ﬂow
137

Water pump
Fan

Motor

Working area

(a)

Working area
IPMC sensor

(b)
Figure 6.7: Loligo Systems swim tunnel: (a) Side view; (b) top view.
velocity using a 515/2536 Rotor-X paddle wheel ﬂow sensor. We mounted the ﬂow sensor
to the swim tunnel as shown in Fig. 6.8. The knob readings of 7 to 30 were used, and the
138

ﬂow was allowed to stabilize at the new ﬂow speed for a minimum of 3 minutes between the
experiments. The calibration curve is shown in Fig. 6.8(c), where we can clearly see an aﬃne
relationship between the sensor output (cm/s) and the knob reading.

Flow sensor

Flow sensor

(a)

(b)
Experiment
Linear fit

Velocity (cm/s)

40
30
20
10
0
0

10

20

Fan speed (r.p.m)

30

(c)
Figure 6.8: Rotor-X ﬂow sensor for swim tunnel calibration: (a) Side view; (b) top view; (c)
calibration curve

6.3.2

Diﬀerentiation of Steady and Unsteady Flow

We ﬁrst mounted an IPMC sensor measuring (long = 3.5 cm × wide = 4 mm × thick =
0.25 mm) and used a Cannon digital camera, to visualize each ﬂow regime experimentally.
139

In the experimental setup, the IPMC sensor was mounted in such a way that it would bend
in the direction parallel to the x-axis. To generate the vortex street, a hollow cylindrical
obstacle with a diameter 6 cm was placed upstream the sensor. The IPMC sensor was
attached, ﬁrst, 4 cm behind the obstacle. Second, the sensor was aligned with the obstacle
edge, which made the cylinder 5 cm upstream from it.
The authors of [5] have reported that, the ﬂow behind obstacle regime can be divided
into three main regions, according to the velocity conditions present in each: (i) base suction
region, (ii) vortex formation region, and (iii) vortex street. First, the average ﬂow frequency
of region (ii) bigger than that region (iii). Second, the suction region (region (i)) has average
ﬂow frequency smaller than those in regions (ii) and (iii). The overall goal is to better
understand the sensing environment from the local perspective of a situated sensor. Three
sensing tasks will be investigated here, diﬀerentiating between steady and unsteady ﬂows,
diﬀerentiating between suction and shed K´rm´n vortex region, and ﬁnally, diﬀerentiating
a a
between big and small size obstacles.
Fig. 6.9 shows the image of an IPMC beam in a uniform ﬂow (in the absence of the
obstacle). We can see that the beam is deﬂected with small-amplitude vibration.
Fig. 6.10 shows that the IPMC sensor is sucked toward the cylinder (predominantly
negative ﬂow) when the sensor is placed right after the cylinder. Then, when the obstacle is
shifted to the right, the sensor bends diﬀerently (predominantly positive ﬂow), as shown in
Fig. 6.11.
This global analysis is useful to reveal that , ﬁrst, the diﬀerences between the uniform
ﬂow and the K´rm´n vortex street can be captured by the sensor beam. Second, it is clear
a a
from Figs. 6.10 and 6.11 that, visually, the IPMC sensor would be able to diﬀerentiate the
140

(a)

(b)

(c)
Figure 6.9: IPMC sensor in a uniform ﬂow in the absence of the obstacle: (a) Sensor at to ;
(b) sensor at t1 , (b) sensor at t2 .

(a)

(b)

Figure 6.10: IPMC sensor right behind obstacle (predominantly negative ﬂow): (a) Illustration; (b) typical sensor conﬁguration.
suction region from the vortex region.
We further examined the IPMC sensor outputs under diﬀerent experimental conditions.
Fig. 6.12 shows a typical sensor response under various ﬂow conditions with a nominal ﬂow
141

(a)

(b)

Figure 6.11: IPMC sensor behind obstacle but shifted from center (predominantly positive
ﬂow): (a) Illustration; (b) typical sensor conﬁguration.
velocity of 10.6 cm/s.

It is clear from Fig. 6.12 that, the water movement behind (at the

center or shifted to the right) a cylinder shows much larger ﬂuctuations than in the free
stream (running water). The frequency of the suction region is diﬀerent (smaller) from the
frequency of the K´rm´n vortex region, as expected.
a a

6.3.3

Diﬀerentiation of Obstacle Size

We tested whether one IPMC sensor could potentially detect and recognize an upstream
cylindrical obstacle. For this purpose, we placed an IPMC sensor ﬁrst behind a cylinder
with larger diameter (6 cm) and then replaced the cylinder with another one that has a
smaller diameter (4 cm), as shown in Fig. 6.13(b). The center coordinates of both cylindrical
obstacles had the same distance away from the IPMC sensor, around 8 cm. The nominal
ﬂow velocity was chosen to be 10.6 cm/s in these experiments.
We can observe from the experiments that, the sensor vibrated with diﬀerent frequencies
and amplitudes when the obstacle size was diﬀerent. Typical IPMC sensor signals for two
cases and the corresponding frequency spectra are shown in Fig. 6.14.
From Fig.6.14, we can see that the vortex shedding frequency does vary with the obstacle
142

size, and the good agreement between the measured and theoretical (Eq. (6.11)) values of
fsh is particularly promising.

6.3.4

Diﬀerentiation of Flow Velocity

We tested whether an IPMC sensor could be used to estimate the free-stream velocity based
on the measurement of vortex shedding frequency. For this purpose, we placed an IPMC
sensor behind a upstream cylinder with a diameter of 6 cm, where the obstacle-sensorseparation distance was 5 cm. We used two nominal ﬂow velocities, 7.5 cm/s and 10.6 cm/s.
Typical IPMC signals are the corresponding frequency spectra are shown in Fig. 6.15.
It can be observed that, when the ﬂow velocity was increased, the vortex-shedding frequency fsh increased as well. The observed trend was consistent with that predicted by
the model (6.11); however, the discrepancy between the measured and predicted values was
relatively large. Quantitative understanding of the relationship between free-stream velocity
and the vortex-shedding frequency needs to be pursued further.

6.3.5

Two IPMC Sensors

We also conducted preliminary experiments on a pair of IPMC sensors, separated out with
a ﬁxed distance (3 cm), to examine the potential of an IPMC array in vortex sensing. In
particular, we tested the capability of the sensors in detecting object movement in air. To
stimulate the two IPMC sensors, a human ﬁnger was moved alongside the sensors at a
distance of 1.5 cm. Fig. 6.16 shows the sensor responses when the ﬁnger was moved in
one direction and then in the opposite direction. It can be observed that there was strong
correlation between the responses of the two sensors, based on which one can compute the
143

time delay between the two responses and consequently the object moving speed.
We conducted further experiments using the two-sensor array, to test whether it could
provide more hydrodynamic information than a single sensor. From those experiments, we
found that: (i) The two sensors were sucked toward the cylinder if both of them were in the
suction region (negatively dominant ﬂow); (ii) they were both bent away from the obstacle
if both of them were in the shed K´rm´n vortex street region (positively dominant ﬂow);
a a
(iii) sensor #1 was sucked toward the obstacle and sensor #2 repelled from the obstacle
when they were in the suction and vortex regions, respectively. Fig. 6.17 shows the case (iii),
where the leading sensor was separated from the obstacle by 4 cm and the nominal velocity
was 10.6 cm/s.
Fig. 6.18 shows typical IPMC signals for the latter scenario, where we can see that the
vortex-shedding frequency in the vortex region is higher than that in the suction region.
These experiments indicate that an IPMC sensor array could potentially be used to map
diﬀerent ﬂow regions in a vortex street.

6.4

Chapter Summary

In this chapter we formulated an estimation problem for K´rm´n vortex street formed by
a a
an obstacle in a ﬂow. While we did not provide a full solution to this problem, a number
of experiments were conducted in a water tunnel with a single or a pair of IPMC sensors,
to shed light on a potential estimation approach. Speciﬁcally, we found that one could use
the vortex-shedding frequency fsh to infer the information about the ﬂow region, obstacle
size, or the nominal ﬂow velocity. The correlation between the signals from multiple sensors
could be used to measure the vortex traveling speed.
144

|Sensor output| (µ A)

Sensor output (µ A)

2
1
0
−1
−2
0

5

10

1

0.5

0
0

2

4

6

Frequency (Hz)

Time (s)

|Sensor output| (µ A)

Sensor output (µ A)

(a)

2
1
0
−1
−2
0

5

10

Time (s)

1

0.5

0
0

5

Frequency (Hz)

3

|Sensor output| (µ A)

Sensor output (µ A)

(b)

2
1
0
−1
−2
0

5

10

Time (s)

1

0.5

0
0

5

Frequency (Hz)
(c)

Figure 6.12: Typical IPMC sensor signals and frequency spectrum of the response under different ﬂow conditions with nominal ﬂow velocity of 10.6 cm/s: (a) Running water (free ﬂow);
(b) predominantly negative ﬂow (suction region), frequency = 0.61 Hz; (c) predominantly
positive ﬂow (vortex region), frequency = 0.67 Hz.
145

(a)

(b)
Figure 6.13: IPMC sensor was attached behind of the: (a) 6 cm diameter cylindrical obstacle;
(b) 4 cm diameter cylindrical obstacle.

146

1
|Sensor output| (µ A)

Sensor output (µ A)

4
3
2
1
0
−1
−2
0

5
Time (s)

10

0.8
0.6
0.4
0.2
0
0

2
4
Frequency (Hz)

6

2
4
Frequency (Hz)

6

(a)

1
|Sensor output| (µ A)

Sensor output (µ A)

3
2
1
0
−1
−2
0

5
Time (s)

10

0.8
0.6
0.4
0.2
0
0

(b)
Figure 6.14: Typical IPMC sensor signals for diﬀerent-size cylindrical obstacles with diameter: (a) 6 cm, theoretical / measured fsh : 0.56/0.6 Hz; (b) 4 cm, theoretical / measured
fsh : 0.73/0.79 Hz.

147

|Sensor output| (µ A)

Sensor output (µ A)

4
2
0
−2
0

5

1
0.8
0.6
0.4
0.2
0
0

10

Time (s)

2

4

6

Frequency (Hz)

3

|Sensor output| (µ A)

Sensor output (µ A)

(a)

2
1
0
−1
−2
0

5

10

Time (s)

1

0.5

0
0

5

Frequency (Hz)
(b)

Figure 6.15: Typical IPMC sensor signals under 6 cm upstream cylinder, the nominal ﬂow
velocity of: (a) 7.5 cm/s, Theoretical / Measured fsh :0.4/0.61 Hz;(b) 10.6 cm/s, Theoretical
/ Measured fsh : 0.56/0.67 Hz.

148

Sensor output (µ A)

1
Sensor #1
Sensor #2
0.5

0

−0.5
0

0.5

1

1.5

2

Time (s)
(a)

Sensor output (µ A)

1
Sensor #1
Sensor #2

0.5
0
−0.5
−1
0

0.5

1

1.5

2

Time (s)
(b)
Figure 6.16: Typical IPMC sensor signals to the human ﬁngers that moved alongside the
two IPMC sensors with velocity of about 25 cm/s, direction of the movement: (a) Toward
positive x-axis; (b) toward negative x-axis.
149

Obstacle

3 cm
4 cm
Sensors

(a)

(b)
Figure 6.17: Two IPMC sensors were attached behind a 6 cm cylinder: (a) Setup at time
to ; (b) sensor #1 sucked toward the obstacle, while sensor#2 repelled from the obstacle at
time t1 .

150

|Sensor output| (µ A)

Sensor output (µ A)

2
0
−2
−4
0

5

10

1

0.5

0
0

Time (s)

2

4

6

Frequency (Hz)

4

|Sensor output| (µ A)

Sensor output (µ A)

(a)

2
0
−2
−4
0

5

10

Time (s)

1

0.5

0
0

2

4

6

Frequency (Hz)
(b)

Figure 6.18: Signals from the two IPMC sensors and their frequency spectra: (a) Sensor #1,
dominant frequency = 0.55 Hz; (b) sensor #2, dominant frequency = 0.73 Hz.

151

Chapter 7
Conclusion
The ﬁrst contribution of this dissertation was a new approach to the realization of artiﬁcial
lateral lines for underwater robots and vehicles that uses IPMC materials as sensing elements.
For example, we showed that, with a neural network-based processing scheme, a dipole
source can be successfully localized when its vibration amplitude and orientation are known.
Experimental results showed that, with relatively few sensor elements, the IPMC-based
lateral line was able to localize sources at least 1-2 BLs away, and the localization accuracy
at source-sensor separation of 1 BL was comparable to the resolution of manually placing
the source (1-2 mm).
Second, we proposed and examined in detail two iterative schemes for dipole source localization based on measurements from an artiﬁcial lateral line. The schemes, referred to
as Gauss-Newton (GN) algorithm and Newton-Raphson (NR) algorithm, respectively, were
derived with a nonlinear estimation perspective to simultaneously estimate the dipole location, vibration amplitude, and vibration orientation. We performed CRB analysis based on
the analytical ﬂow model, and demonstrated its use in the design of the lateral line system
152

via the example of optimizing intra-sensor spacing. The CRB analysis was also used later
to conﬁrm the near-optimality of the proposed algorithms. We conducted both simulation
and physical experimentation to validate the proposed algorithms. The good consistency between the simulation results and the experimental results does not only suggest the physical
relevance of our simulation setup, but also supports the validity of the analytical models and
algorithms for the physical system. With a body length (BL) of 10 cm for the lateral line
system, simulation and experimental results showed that both the GN and NR algorithms
can successfully localize the dipole source and estimate its vibration amplitude and orientation within 1 BL, with the maximum localization error less than 5% of the BL. We also
found that, while the GN algorithm did consistently better than the NR algorithm, the overall performance and computational complexity of the two algorithms are very comparable.
We also found that the proposed algorithms showed clear advantage over the beamforming
algorithm in terms of accuracy and computational eﬃciency.
Third, we investigated an analytical model-based estimation approach to the localization
of a moving dipole source using an artiﬁcial lateral line consisting of an array of ﬂow sensors.
The moving speed, vibration amplitude, and location of the moving source were assumed to
be unknown, and we proposed a Gauss-Newton algorithm to solve the nonlinear estimation
problem based on the evolving signal magnitudes obtained with SDFT. The eﬀectiveness of
the approach was illustrated with simulation and experimental results.
Fourth, a novel two-stage extended Kalman ﬁlter was proposed to track and estimate
a moving but non-vibrating object in water. Simulation and experimental results have
demonstrated that, with a few IPMC sensors, the ALL was able to identify and estimate the
moving object speed, location, size, and shape parameter, respectively.
153

Fifth, we conducted experimental exploration of using IPMC sensors for detecting hydrodynamic parameters associated with a K´rm´n vortex street. In particular, our work
a a
showed that it was promising to use IPMC lateral line for the estimation of obstacle size,
free-stream ﬂow velocity, distance to the obstacle, etc.
Although this work has demonstrated the promise of IPMC sensors in underwater ﬂow
sensing, we note that the behavior of IPMCs depends on environmental conditions and
could vary over time. While the continuing advances in IPMC fabrication may alleviate such
problems, appropriate compensation schemes (see, e.g., [35]) and periodic sensor calibrations
are expected to be essential for practical IPMC lateral line systems.

154

Chapter 8
Future Work
The work reported in this dissertation can be extended in several directions. First, in this
work we largely ignored the eﬀect of sensor beams on the ﬂow ﬁeld. On the sensor modeling
side, we primarily took the short-circuit current output of the sensor as proportional to the
ﬂow velocity. It will be of interest to extend the modeling work on the sensor-ﬂow interaction,
which can potentially lead to improved sensing performance for the artiﬁcal lateral line.
Second, while we had demonstrated the eﬀectiveness of a two-stage extended Kalman
ﬁltering scheme in object tracking and estimation, the convergence and stability properties
of such schemes were not examined in this dissertation. The study of these issues, possibly
in a more general n-stage setting (Fig. 8.1), would be of interest.
Third, while we formulated an estimation problem for the case of vortex street sensing,
we did not provide a systematic solution to this problem. It would be interesting to develop
a complete, analytical approach for the estimation of obstacle location, size, ﬂow speed, and
ﬂow direction, and validate the approach experimentally.
Fourth, comparing to the size (less than 1 mm long) of neuromasts in a biological lateral
155

Figure 8.1: Illustration of the n-stages nonlinear ﬁltering scheme.

line, the IPMC sensors used in this dissertation are relatively big. In addition, a biological
lateral line may have many more (from tens to over a thousand [20]) neuromasts. There
are many challenges in realizing an IPMC-based artiﬁcial lateral line system with sensor
size and number comparable to those of a live ﬁsh. For example, while the processes for
micro fabrication of IPMCs have been reported by a number of groups [18, 101, 32], existing
results deal mostly with planar processes, and new polymer MEMS techniques have to be
developed to create micro IPMC devices standing on a substrate as in [54]. On the system
integration side, the processing of signals from tens to hundreds of IPMC sensors is far from
trivial. It will be of interest to extend the signal processing approaches developed in this
dissertation to accommodate a much larger number of sensors. To address this challenge, it
156

will be critical to seek inspiration from how such processing is achieved in biology [21].
Finally, artiﬁcial lateral line system can be mounted on a robot, such as a robotic ﬁsh,
for the investigation of processing schemes that decouple external signals from self-motioninduced ﬂow signals. Such an extension would be particularly relevant for practical applications of artiﬁcial lateral lines.

157

APPENDICES

158

Appendix A

Artiﬁcial Lateral Line Design: CRB

With θ = [xs , ys ], the ﬁrst and second derivatives are given as:
∂fi (θ)
=
∂xs

k1
ri 5
5k1
ri

7

[−4(xi − xs )α1 − 3(yi − ys )α2 ] +
(2(xi − xs )3 − (xi − xs )(yi − ys )2 )α1 + (3(xi − xs )2 (yi − ys ))α2 ,

5k1
∂ 2 fi (θ)
4k1α1
+
(−10(xi − xs )2 + (yi − ys )2 )α1 − (9(xi − xs )(yi − ys ))α2
=
5
2
(∂xs )
ri
ri 7
35k1
(2(xi − xs )4 − (xi − xs )2 (yi − ys )2 )α1 + (3(xi − xs )3 (yi − ys ))α2 ,
9
ri

∂fi (θ)
=
∂ys

k1
ri 5
5k1

[2(yi − ys )α1 − 3(xi − xs )α2 ] +

(2(xi − xs )2 (yi − ys ) − (yi − ys )3 )α1 + (3(xi − xs )(yi − ys )2 )α2 ,
ri 7
159

−2k1α1
∂ 2 fi (θ)
5k1
=
+
(−2(xi − xs )2 + 5(yi − ys )2 )α1 + (−9(xi − xs )(yi − ys ))α2
2
5
7
(∂ys )
ri
ri
35k1
(2(xi − xs )2 (yi − ys )2 − (yi − ys )4 )α1 + (3(xi − xs )(yi − ys )3 )α2 ,
9
ri

∂ 2 fi (θ)
∂ 2 fi (θ)
=
∂ys ∂xs
∂xs ∂ys
5k1
3k1α2
+
(−2(xi − xs )(yi − ys ))α1 − 3((xi − xs )2 + (yi − ys )2 )α2
=
5
ri
ri 7
35k1
(2(xi − xs )3 (yi − ys ) − (xi − xs )(yi − ys )3 )α1
9
ri
+(3(xi − xs )2 (yi − ys )2 )α2 ,
where k1 = a3 πf .

160

(A.1)

Appendix B
Appendix for Chapter5

B.1

Derivation of CRB

Based on the measurement equation (5.11), the likelihood function is given by
2
N −1
i=0 (Mi −fi (X))

−1

1
2
p(M ; X) =
e 2πσm
2
2πσm

.

(B.1)

Diﬀerentiating the log of the likelihood function twice, we get:∂2
ln (p(M ; X)) =
∂ 2X

1
2
σm

N −1
i=0

(Mi − fi (X))

∂ 2 fi (X)
∂fi (X)
−
2X
∂X
∂

The Fisher information matrix is deﬁned by −E
1
F = 2
σm

N −1
i=0

∂fi (X)
∂X
161

∂fi (X) T
∂X

∂ 2 (ln p(M ;X))
, which is given by
(∂X)2

∂fi (X) T
∂X

.

(B.2)

where E[·] is the expectation taken over all possible measurement M . the ith diagonal of
−1
the inverse of F , which is denoted by Fi,i , provides a lower bound on the variance of the
−1
ˆ
ˆ
ith element of the state estimate X. Generally, V ar(X) ≥ Fi,i . In our case we have four

elements and F is just a 4 × 4 matrix. The Fisher information matrix can be expressed as
T
F = [F11 , F12 ; F12 , F22 ], where the sub-matrix blocks F11 , F12 , and F22 are 2 × 2 matrices,



N
i=1

∂fi (·) 2
∂xs

N
i=1



∂fi (·)
∂ys 
,
∂fi (·) 2
ys

∂fi (·)
∂xs

1 
F11 = 2 
σm

N
i=1

∂fi (·)
ys

∂fi (·)
xs



N
i=1

∂fi (·)
∂xs

∂fi (·)
∂R

N
i=1

∂fi (·)
∂xs

N
i=1

∂fi (·)
∂ys

∂fi (·)
∂R

N
i=1

∂fi (·)
∂ys

N
i=1

∂fi (·)
∂R

1 
F12 = 2 
σm


1 
F22 = 2 
σm

N
i=1
N
i=1

N
i=1

∂fi (·) 2
∂R

∂fi (·)
∂λ1

∂fi (·)
∂R

N
i=1



∂fi (·)
∂λ1 

∂fi (·)
∂λ1

,



∂fi (·)
∂λ1 
,
∂fi (·) 2
∂λ1

where

∂f (·)
= vx R 2
∂xs
−

2(x − xs )

(x − xs )2 − (y − ys )2 − λ1

2

+ 4(x − xs )2 (y − ys )2

4(x − xs ) (x − xs )2 − (y − ys )2 − λ1

(x − xs )2 + (y − ys )2 − λ1

2
2
(x − xs )2 − (y − ys )2 − λ1 + 4(x − xs )2 (y − ys )2

162







,

∂f (·)
= vx R 2
∂ys
+

−2(y − ys )

(x − xs )2 − (y − ys )2 − λ1

2

+ 4(x − xs )2 (y − ys )2




−(x − xs )2 − (y − ys )2 − λ1 
,
2

2

(x − xs )2 − (y − ys )2 − λ1 + 4(x − xs )2 (y − ys )2

4(y − ys ) (x − xs )2 − (y − ys )2 − λ1

∂f (·)
= −2vx R
∂R

∂f (·)
= vx R 2
∂λ1
−

(x − xs )2 − (y − ys )2 − λ1

(x − xs )2 − (y − ys )2 − λ1

2

+ 4(x − xs )2 (y − ys )2

1
2

+ 4(x − xs )2 (y − ys )2


2

2 (x − xs )2 − (y − ys )2 − λ1
,
2

2 − (y − y )2 − λ 2 + 4(x − x )2 (y − y )2 
(x − xs )
s
s
s
1
(x − xs )2 − (y − ys )2 − λ1

and the Cramer-Rao bound is deﬁned as F −1 , where F −1 is given by



F −1 = 

−1
−JF12 F22

J −1

−1
−1
F22 F21 J −1 F22 (I + F21 J −1 F12 )

−1
where J = F11 − F12 F22 F21 . J is a 2 × 2 matrix given by





J1,1 J1,2 
J =
,
J2,1 J2,2

163




,

,

−1
−1
The CRBs for the estimates of xs and ys are given by J1,1 = |J|−1 ·J2,2 and J2,2 = |J|−1 ·J1,1 ,

respectively. Speciﬁcally,

−1
−1
CRB(xs ) = J1,1 , CRB(ys ) = J2,2 .

B.2

(B.3)

More Results for One-stage Extended Kalman Filtering Scheme

In this appendix we present similar results as in Section 5.3 but with a one-stage extended
Kalman ﬁltering scheme. First, we have tested the impact of the process noise on the onestage scheme as shown in Fig. B.1.
Table B.1 summarizes the other estimated parameters. It is clear from the simulation
results of the one-stage scheme that, the model process noise has slight impact on the estimation results.
Table B.1: One-stage scheme: Estimated velocity, size parameter, and shape parameter
under diﬀerent process noise variances.

Parameter
Velocity (cm/s)
Size R (cm)
Shape λ1 (cm)

Actual
6
2.5
1.5

Estimated (σ 2 = 0.1)
6.33
2.45
0.63

(σ 2 = 0.3)
6.53
2.45
1.66

(σ 2 = 0.5)
6.56
2.12
0.71

Fig. B.2 shows the impact of the measurement noises on the estimation localization results
when the process noise is held at 0.3. Table B.2 lists the other estimated parameters. It can
be seen that, overall, the measurement noise variance has a large impact on the estimation
performance; in general, larger measurement noise variance leads to much longer transient
164

Actual track
Estimated track (σ2 = 0.1)
Estimated track (σ2 = 0.3)
Estimated track (σ2 = 0.5)

y (cm)

8
6
4
2
0
0

2

4

6

8

10

x (cm)
(a)
2

σ = 0.1
σ2 = 0.3
σ2 = 0.5

Localization error (cm)

1
0.8
0.6
0.4
0.2
0
0

2

4

6

8

10

x (cm)
(b)
Figure B.1: Simulation results on tracking the ellipsoidal cylinder for one-stage scheme: (a)
The estimation of the moving object location for diﬀerent variance values for the process
noises; (b) localization error at each x-coordinate of the object. The measurement noise
variance is held at 0.01 cm2 /s2 .

165

Actual track
2

Estimated track (σm = 0.01)
2

Estimated track (σm = 0.03)
2
m

Estimated track (σ = 0.05)

y (cm)

8
6
4
2
0
0

2

4

6

8

10

x (cm)
(a)
2

σm = 0.01
2

σm = 0.03

Localization error (cm)

1.5

2

σm = 0.05

1

0.5

0
0

2

4

6

8

10

x (cm)
(b)
Figure B.2: Simulation results on tracking the ellipsoidal cylinder for one-stage scheme: (a)
The estimation of the moving object location for diﬀerent variance values for the measurement noises; (b) localization error at each x-coordinate of the object. The process noise
variance is held at 0.3 (in its respective unit for each state component).

166

and larger estimation error.
Table B.2: One-stage scheme: Estimated velocity, size parameter, and shape parameter
under diﬀerent measurement noise variances.
Parameter
Velocity( cm/s)
Size R (cm)
Shape λ1 (cm)

Actual
6
2.5
1.5

2
Estimated (σm = 0.01)
6.45
2.33
0.85

167

2
(σm = 0.03)
6.63
2.15
0.58

2
(σm = 0.05)
7.23
2.64
0.42

BIBLIOGRAPHY

168

BIBLIOGRAPHY
[1] A. T. Abdulsadda and X. Tan. Underwater source localization using an IPMC-based
artiﬁcial lateral line. In Proceedings of the 2011 IEEE International Conference on
Robotics and Automation, pages 447–452, Shanghai, China, 2011.
[2] A. T. Abdulsadda and X. Tan. An artiﬁcial lateral line system using IPMC sensor arrays. International Journal of Smart and Nano Materials, 2012. to appear,
DOI:10.1080/19475411.2011.650233.
[3] A. T. Abdulsadda, F. Zhang, and X. Tan. Localization of source with unknown amplitude using IPMC senosr arrays. In Y. Bar-Cohen, editor, Proceedings of Electroactive
Polymer Actuators and Devices (EAPAD) XIII, volume 7976 of Proceedings of SPIE,
page 797627, San Diego, CA, 2011.
[4] O. Akanyeti, R. Venturelli, F. Visentin, L. Chambers, W. Megill, and P. Fiorini. What
information do k´rm´n streets oﬀer to ﬂow sensing? Bioinsp. Biomim., 6:036001,
a a
2011.
[5] Otar Akanyeti, Roberto Venturelli, Francesco Visentin, Lily Chambers William Megill,
and Paolo Fiorini. What information do k´rm´n streets oﬀer to ﬂow sensing? Bioina a
spiration and Biomimetics, 6:036001:12, 2011.
[6] B. Akle, M. Bennett, and D. Leo. High-strain ionomeric-ionic liquid electroactive
actuators. Sensors and Actuators A, 126:173–181, 2006.
[7] M. Aureli, V. Kopman, and M. Porﬁri. Free-locomotion of underwater vehicles actuated by ionic polymer metal composites. IEEE/ASME Transactions on Mechatronics,
15(4):603–614, 2010.
[8] M. Aureli, C. Prince, M. Porﬁri, and S. D. Peterson. Energy harvesting from base
excitation of ionic polymer metal composites in ﬂuid environments. Smart Materials
and Structures, 19:015003, 2010.
169

[9] H. Bleckmann. Peripheral and central processing of lateral line information. J. Comp.
Physiol. A, 194:145–158, 2008.
[10] C. Bonomo, L. Fortuna, P. Giannone, S. Graziani, and S. Strazzeri. A model for ionic
polymer metal composites as sensors. Smart Materials and Structures, 15:749–758,
2006.
[11] R. Bouﬀanais, G. D. Weymouth, and D. K. P. Yue. Hydrodynamic object recognition
using pressure sensing. Proc. R. Soc. A, 467:19–38, 2011.
[12] Paola Brunetto, Luigi Fortuna, Salvatore Graziani, and Salvatore Strazzeri. A model
of ionic polymer-metal composite actuators in underwater operations. Smart Materials
and Structures, 17(2):025029:1–12, 2008.
[13] X. Chen, G. Zhu, X. Yang, D. L. S. Hung, and X. Tan. Model-based estimation
of ﬂow characteristics using an ionic polymermetal composite beam. IEEE/ASME
Transactions on Mechatronics, 2012. to appear, DOI: 10.1109/TMECH.2012.2194300.
[14] Z. Chen, D. Hedgepeth, and X. Tan. A nonlinear, control-oriented model for ionic
polymer-metal composite actuators. Smart Materials and Structures, 18:055008:1–9,
2009.
[15] Z. Chen, S. Shatara, and X. Tan. Modeling of biomimetic robotic ﬁsh propelled by an
ionic polymer-metal composite caudal ﬁn. IEEE/ASME Transactions on Mechatronics,
15(3):448–459, 2010.
[16] Z. Chen, Y. Shen, N. Xi, and X. Tan. Integrated sensing for ionic polymer-metal composite actuators using PVDF thin ﬁlms. Smart Materials and Structures, 16(2):S262–
S271, 2007.
[17] Z. Chen and X. Tan. A control-oriented and physics-based model for ionic polymermetal composite actuators. IEEE/ASME Transactions on Mechatronics, 13(5):519–
529, 2008.
[18] Z. Chen and X. Tan. Monolithic fabrication of ionic polymer-metal composite actuators
capable of complex deformation. Sensors and Actuators A: Physical, 157:246–257,
2010.
[19] Z. Chen, X. Tan, A. Will, and C. Ziel. A dynamic model for ionic polymer-metal
composite sensors. Smart Materials and Structures, 16:1477–1488, 2007.
170

[20] S. Coombs. Smart skins: Information processing by lateral line ﬂow sensors. Autonomous Robots, 11:255–261, 2001.
[21] S. Coombs and C. B. Braun. Information processing by the lateral line system. In
S. P. Collin and N. J. Marshall, editors, Sensory Processing in Aquatic Environments,
pages 122–138, New York, 2003. Springer-Verlag.
[22] S. Coombs and R. A. Conley. Dipole source localization by mottled sculpin. I. Approach
strategies. J. Comp. Physiol. A, 180:387–399, 1997.
[23] S. Coombs and R. A. Conley. Dipole source localization by mottled sculpin. II. the
role of lateral line excitation patterns. J. Comp. Physiol. A, 180:401–415, 1997.
[24] S. Coombs, M. Hastings, and J. Finneran. Modeling and measuring lateral line excitation patterns to changing dipole source locations. J. Comp. Physiol. A, 178:359–371,
1996.
[25] T. M. Cover and J. A. Thomas. Elements of Information Theory. John Wiley & Sons,
Inc., New York, 1991.
[26] B. Curcic-Blake and S. M. van Netten. Source location encoding in the ﬁsh lateral line
canal. Journal of Experimental Biology, 209:1548–1559, 2006.
[27] A. M. K. Dagamseh, T. S. J. Lammerink, M. L. Kolster, C. M. Bruinink, R. J.
Wiegerink, and G. J. M. Krijnen. Dipole-source localization using biomimetic ﬂowsensor arrays positioned as lateral-line system. Sensors and Actuators A, 162(2):355–
360, 2010.
[28] P. de Gennes, K. Okumura, M. Shahinpoor, and K. Kim. Mechanoelectric eﬀects in
ionic gels. Europhys. Lett., 50(4):513–518, 2000.
[29] M. Dijkstra, J. J. van Baar, R. J. Wiegerink, T. S. J. Lammerink, J. H. de Boer, and
G. J. M. Krijnen. Artiﬁcial sensory hairs based on the ﬂow sensitive receptor hairs of
crickets. J. Micromech. Microeng., 15:S132–S138, 2005.
[30] A. F. Emery and A. V. Nenarokomov. Optimal experiment design. Measurement
Science and Technology, 9(6):864–876, 1998.
[31] Z. Fan, J. Chen, J. Zou, D. Bullen, C. Liu, and F. Delcomyn. Design and fabrication
of artiﬁcial lateral line ﬂow sensors. Journal of Micromechanics and Microengineering,
12:655–661, 2002.
171

[32] G.-H. Feng and R.-H. Chen. Fabrication and characterization of arbitrary shaped
µ IPMC transducers for accurately controlled biomedical applications. Sensors and
Actuators:A physical, 143(1):34–40, 2008.
[33] V. I. Fernandez, A. Maertens, F. M. Yaul, J. Dahl, J. H. Lang, and M. S. Triantafyllou.
Lateral-line-inspired sensor arrays for navigation and object identiﬁcation. Marine
Technology Socity Journal, 45(4):130–146, 2011.
[34] J.-M. P. Franosch, A. B. Sichert, M. D. Suttner, and J. L. van Hemmen. Estimating
position and velocity of a submerged moving object by the clawed frog xenopus and
by ﬁsh: A cybernetic approach. Biol. Cybern., 93:231–238, 2005.
[35] T. Ganley, D. L.S. Hung, G. Zhu, and X. Tan. Modeling and inverse compensation of
temperature-dependent ionic polymer-metal composite sensor dynamics. IEEE/ASME
Transactions on Mechatronics, 16(1):80–89, 2011.
[36] J. M. Gardiner and J. Atema. Sharks need the lateral line to locate odor sources:
Rheotaxis and eddy chemotaxis. Journal of Experimental Biology, 210:1925–1934,
2007.
[37] Stokes G.G. On the eﬀect if the internal friction of ﬂuids on the motion of pendulums.
Trans. Camb. Philos Soc., 1201:9:8–106, 1951.
[38] Batchelor G.K. An introduction to ﬂuid dynamics. Cambridge University Press, New
York, 1967.
[39] J. Goulet, J. Engelmann, B. P. Chagnaud, J.-M. P. Franosch, M. D. Suttner, and J. L.
van Hemmen. Object localization through the lateral line system of ﬁsh: theory and
experiment. J. Comp. Physiol. A, 194:1–17, 2008.
[40] E.-S. Hassan. Mathematical analysis of the stimulus for the lateral line organ. Biol.
Cybern., 52:23–36, 1985.
[41] E.-S. Hassan. Mathematical description of the stimuli to the lateral line system of ﬁsh
derived from a three-dimensional ﬂow ﬁeld analysis: I. The cases of moving in open
water and of gliding towards a plane surface. Biol. Cybern., 66:443–452, 1992.
[42] E.-S. Hassan. Mathematical description of the stimuli to the lateral line system of ﬁsh
derived from a three-dimensional ﬂow ﬁeld analysis: II. The case of gliding alongside
or above a plane surface. Biol. Cybern., 66:453–461, 1992.
172

[43] E.-S. Hassan. Mathematical description of the stimuli to the lateral line system of
ﬁsh derived from a three-dimensional ﬂow ﬁeld analysis: III. The case of an oscillating
sphere near the ﬁsh. Biol. Cybern., 69:525–538, 1993.
[44] E. Jacobsen and R. Lyons. The sliding dft. IEEE Signal Process. Mag., 20(2):74–80,
2003.
[45] Gere J.M. Mechanics of materials. Nelson Thornes Ltd, Cheltenham, 2001.
[46] A. J. Kalmijn. Hydrodynamic and acoustic ﬁeld detection. In J. Atema, R. R. Fay,
A. N. Popper, and W. N. Tavolga, editors, Sensory Biology of Aquatic Animals, pages
83–130, New York, 1988. Springer-Verlag.
[47] S. M. Kay. Fundamentals of Statistical Signal Processing: Estimation Theory. Prentice
Hall, Upper Saddle River, New Jersey, 1993.
[48] K. J. Kim, D. Pugal, and K. K. Leang. A twistable ionic polymer-metal composite
artiﬁcial muscle for marine applications. Marine Technology Society Journal, 45(4):83–
98, 2011.
[49] K.J. Kim and M. Shahinpoor. Ionic polymer-metal composites: II. manufacturing
techniques. Smart Materials and Structures, 12:65–79, 2003.
[50] A. Klein and H. Bleckmann. Determination of object position, vortex shedding frequency and ﬂow velocity using artiﬁcial lateral line canals. Beilstein Journal of Nanotechnology, 2:276–283, 2011.
[51] A. Kong, J. S. Liu, and W. H. Wong. Sequential imputations and bayesian missing
data problem. Journal of the American Statistical Association, 89(425):278–288, 1994.
[52] A. B. A. Kroese and N. A. M. Schellart. Velocity- and acceleration-sensitive units in
the trunk lateral line of the trout. Journal of Neurophysiology, 68(6):2212–2221, 1992.
[53] H. Lamb. Hydrodynamics. Cambridge University Press, Cambridge, 1932.
[54] H. Lei, W. Li, and X. Tan. Microfabrication of IPMC cilia for bio-inspired ﬂow sensing.
In Y. Bar-Cohen, editor, Electroactive Polymer Actuators and Devices (EAPAD) XIV,
volume 8340 of Proceedings of SPIE, page 83401A, Bellingham, WA, 2012. SPIE.
173

[55] J. C. Liao. Neuromuscular control of trout swimming in a vortex street: Implications for energy economy during the k´rm´n gait. Journal of Experimental Biology,
a a
207:3495–3506, 2004.
[56] J. C. Liao, D. N. Beal, G. V. Lauder, and M. S. Triantafyllou. Fish exploiting vortices
decrease muscle activity. Science, 302:1566–1569, 2003.
[57] James C Liao. A review of ﬁsh swimming mechanics and behaviour in altered ﬂows.
Phil. Trans. Royal Society, 2082:1973–1993, 2007.
[58] C. Lim, H. Lei, and X. Tan. A dynamic, physics-based model for base-excited IPMC
sensors. In Y. Bar-Cohen, editor, Electroactive Polymer Actuators and Devices (EAPAD) XIV, volume 8340 of Proceedings of SPIE, page 83400H, Bellingham, WA, 2012.
SPIE.
[59] C. Liu. Micromachined biomimetic artiﬁcial haircell sensors.
2:S162–S169, 2007.

Bioinsp. Biomim.,

[60] M. E. McConney, N. Chen, D. Lu, H. A. Hu, S. Coombs, C. Liu, and V. V. Tsukruk.
Biologically inspired design of hydrogel-capped hair sensors for enhanced underwater
ﬂow detection. Soft Matter, 5:292–295, 2009.
[61] M. J. McHenry, J. A. Strother, and S. M. van Netten. Mechanical ﬁltering by the
boundary layer and ﬂuidstructure interaction in the superﬁcial neuromast of the ﬁsh
lateral line system. J. Comp. Physiol. A, 194:795–810, 2008.
[62] C. Musso, N. Oudjane, and F. LeGland. Improving regularity particle ﬁlters. SpringerVerlag, New York, 2001.
[63] S. Nemat-Nasser and J. Li. Electromechanical response of ionic polymermetal composites. J. Appl.Phys., 87(7):3321–3331, 2000.
[64] NeuroSolutions. http://www.neurosolutions.com/.
[65] Nam Nguyen, Douglas L. Jones, Yingchen Yang, and Chang Liu. Flow vision for
autonomous underwater vehicles via an artiﬁcial lateral line. EURASIP Journal on
Advances in Signal Processing, 2011:806406:1–11, 2011.
[66] A. V. Oppenheim and R. W. Schafer. Discrete-Time Signal Processing. Pearson Higher
Education, Inc., Upper Saddle River, New Jersey, 3rd edition, 2010.
174

[67] Y. Ozaki, T. Ohyama, T. Yasuda, and I. Shimoyama. An air ﬂow sensor modeled on
wind receptor hairs of insects. In Proceedings of the 13th IEEE International Conference on Micro Electro Mechanical Systems, pages 531–536, Miyazaki, Japan, 2000.
[68] S. Pandya, Y. Yang, D. L. Jones, J. Engel, and C. Liu. Multisensor processing algorithms for underwater dipole localization and tracking using MEMS artiﬁcial lateralline sensors. EURASIP Journal on Applied Signal Processing, 2006, 2006. Article ID
76593, 8 pages, 2006. doi:10.1155/ASP/2006/76593.
[69] T.J. Pitcher, B. L. Partridge, and C. S. Wardle. A blind ﬁsh can school. Science,
194:963–965, 1976.
[70] K. Pohlmann, J. Atema, and T. Breithaupt. The importance of the lateral line in
nocturnal predation of piscivorous catﬁsh. Journal of Experimental Biology, 207:2971–
2978, 2004.
[71] M. Porﬁri. Charge dynamics in ionic polymer metal composites. J. Appl.Phys.,
104(10):104915:1–10, 2008.
[72] D. Pugal, K. J. Kim, and A. Aabloo. Modeling the transduction of IPMC in 3D conﬁgurations. In Behavior and Mechanics of Multifunctional Materials and Composites
2010, volume 7644 of Proceedings of SPIE, page 76441T, San Diego, CA.
[73] D. Pugal, K.Jung, A. Aabloo, and K. J. Kim. Ionic polymer-metal composite mechanoelectrical transduction: Review and perspectives. Polymer International, 59:279–289,
2010.
[74] A. Qualtieri, F. Rizzi, M.T. Todaro, A. Passaseo, R. Cingolani, and M. De Vittorio.
Stress-driven AlN cantilever-based ﬂow sensor for ﬁsh lateral line system. Microelectronic Engineering, 88(8):2376–2378, 2011.
[75] Z. Ren and K. Mohseni. A model of the lateral line of ﬁsh for vortex sensing. Bioinspir.
Biomim., 7:036016, 2012.
[76] Branko Ristic, Sanjeev Arulampalam, and Neil Gordon. Beyond the Kalman Filter.
Artech House, London, 2004.
[77] Panton RL. Incompressible ﬂow. Wiley, New York, 1985.
175

[78] M. M. Sadeghi, R. L. Peterson, and K. Najaﬁ. Micro-hydraulic structure for high
performance bio-mimetic air ﬂow sensor arrays. In Proceedings of the 2011 IEEE
International Electron Devices Meeting, pages 29.4.1–29.4.4, 2011.
[79] Sandhya Samarasinghe. Neural Networks for Applied Sciences and Engineering: From
Fundamentals to Complex Pattern Recognition. Auerbach, Boca Raton, FL, 2007.
[80] S. A. Sarles, P. Pinto, and D. J. Leo. Hair cell sensing with encapsulated interface
bilayers. In Proceedings of Bioinspiration, Biomimetics, and Bioreplication, volume
7975 of Proceedings of SPIE, page 797509, San Diego, CA, 2011.
[81] M. Shahinpoor and K.J. Kim. Ionic polymer-metal composites: I. Fundamentals.
Smart Materials and Structures, 10:819–833, 2001.
[82] S. Shatara and X. Tan. An eﬃcient, time-of-ﬂight-based underwater acoustic ranging
system for small robotic ﬁsh. IEEE Journal of Oceanic Engineering, 35(4):837846,
2010.
[83] A. B. Sichert, R. Bamler, and J. L. van Hemmen. Hydrodynamic object recognition:
When multipoles count. Physical Review Letters, 102:058104, 2009.
[84] A. B. Sichert and J. L. van Hemmen. How stimulus shape aﬀects lateral-line perception:
analytical approach to analyze natural stimuli characteristics. Biol. Cybern., 102:177–
180, 2010.
[85] B.W. Silverman. Density Estimation for Statistics and Data Analysis. Chapman and
Hall, 1986.
[86] K. Takagi, Y. Nakabo, Z. Luo, T. Mukai, M. Yamamura, and Y. Hayakawa. An
analysis of the increase of bending response in IPMC dynamics given uniform input.
In Smart Structures and Materials 2006: Electroactive Polymer Actuators and Devices
(EAPAD), volume 6168 of Proceedings of SPIE, page 616814, San Diego, CA.
[87] K. Takagi, M. Yamamura, Z. Luo, M. Onish, S. Hirano, K. Asaka, and Y. Hayakawa.
Development of a rajiform swimming robot using ionic polymer artiﬁcial muscles. In
Proceedings of the 2006 IEEE/RSJ International Conference on Intelligent Robots and
Systems, pages 1861–1866, Beijing, China, 2006.
[88] Xiaobo Tan. Autonomous robotic ﬁsh as mobile sensor platforms: Challenges and
potential solutions. Marine Technology Society Journal, 45(4):31–40, 2011.
176

[89] Milne Thomson. Theoretical hydrodynamics. Macmillan Company, New York, USA,
5th edition, 1955.
[90] M. S. Triantafyllou, D. K. P. Yue, and G. S. Triantafyllou. Hydrodynamics of ﬁshlike
swimming. Annu. Rev. Fluid Mech., 32:33–53, 2000.
[91] R. Venturelli, O. Akanyeti, F. Visentin, J. Jezov, L. D. Chambers, G. Toming,
J. Brown, M. Kruusmaa, W. M. Megill, and P. Fiorini. Hydrodynamic pressure sensing with an artiﬁcial lateral line in steady and unsteady ﬂows. Bioinspir. Biomim.,
7:036004, 2012.
[92] Y. Wang, C. Lee, and C. Chiang. A MEMS-based air ﬂow sensor with a free-standing
microcantilever structure. Sensors, 7:2389–2401, 2007.
[93] D. Weihs. Hydrodynamics of ﬁsh schooling. Nature, 241:290–291, 1973.
[94] S. P. Windsor, D. Tan, and J. C. Montgomery. Swimming kinematics and hydrodynamic imaging in the blind Mexican cave ﬁsh (astyanax fasciatus). Journal of Experimental Biology, 211:2950–2959, 2008.
[95] T. Y. Wu. Hydromechanics of swimming propulsion. Part 1. Swimming of a twodimensional ﬂexible plate at variable forward speeds in an inviscid ﬂuid. Journal of
Fluid Mechanics, 46:337–355, 1971.
[96] Y. Yang, J. Chen, J. Engel, S. Pandya, N. Chen, C. Tucker, S. Coombs, D. L. Jones, and
C. Liu. Distant touch hydrodynamic imaging with an artiﬁcial lateral line. Proceedings
of the National Academy of Sciences, 103(50):18891–18895, 2006.
[97] Y. Yang, A. Klein, H. Bleckmann, and C. Liu. Artiﬁcial lateral line canal for hydrodynamic detection. Applied Physics Letters, 99:023701, 2011.
[98] Y. Yang, N. Nguyen, N. Chen, M. Lockwood, C. Tucker, H. Hu, H. Bleckmann, C. Liu,
and D. L. Jones. Artiﬁcial lateral line with biomimetic neuromasts to emulate ﬁsh
sensing. Bioinsp. Biomim., 5:016001, 2010.
[99] W. Yim, J. Lee, and K. J. Kim. An artiﬁcial muscle actuator for biomimetic underwater
propulsors. Bioinspiration & Biomimetics, 2:S31–S41, 2007.
[100] U. Zangrilli and L. M. Weiland. Prediction of the ionic polymer transducer sensing of
shear loading. Smart Materials and Structures, 20:094013:1–9, 2011.
177

[101] J. W. L. Zhou, H.-Y. Chan, A. K. H. To, K. W. C. Lai, and W. J. Li. Polymer
mems actuators for underwater micromanipulation. IEEE/ASME Transactions on
Mechatronics, 9(2):334–342, 2004.

178