Aircraft tracking and classification with VHF passive bistatic radar

Abstract

Piotr Ptak
AIRCRAFT TRACKING AND CLASSIFICATION WITH VHF PASSIVE BISTATIC RADAR
Lappeenranta, 2015
118 p.
Acta Universitatis Lappeenrantaensis 647
Diss. Lappeenranta University of Technology
ISBN 978-952-265-815-9, ISBN 978-952-265-816-6 (PDF), ISSN-L 1456-4491, ISSN 1456-4491

Since the times preceding the Second World War the subject of aircraft tracking has been a core in-
terest to both military and non-military aviation. During subsequent years both technology and con-
ﬁguration of the radars allowed the users to deploy it in numerous ﬁelds, such as over-the-horizon
radar, ballistic missile early warning systems or forward scatter fences. The latter one was arranged
in a bistatic conﬁguration. The bistatic radar has continuously re-emerged over the last eighty years
for its intriguing capabilities and challenging conﬁguration and formulation. The bistatic radar ar-
rangement is used as the basis of all the analyzes presented in this work.

The aircraft tracking method of VHF Doppler-only information, developed in the ﬁrst part of this
study, is solely based on Doppler frequency readings in relation to time instances of their appear-
ance. The corresponding inverse problem is solved by utilising a multistatic radar scenario with
two receivers and one transmitter and using their frequency readings as a base for aircraft trajectory
estimation. The quality of the resulting trajectory is then compared with ground-truth information
based on ADS-B data.

The second part of the study deals with the developement of a method for instantaneous Doppler
curve extraction from within a VHF time-frequency representation of the transmitted signal, with
a three receivers and one transmitter conﬁguration, based on a priori knowledge of the probability
density function of the ﬁrst order derivative of the Doppler shift, and on a system of blocks for
identifying, classifying and predicting the Doppler signal. The extraction capabilities of this set-up
are tested with a recorded TV signal and simulated synthetic spectrograms.

Further analyzes are devoted to more comprehensive testing of the capabilities of the extraction
method. Besides testing the method, the classiﬁcation of aircraft is performed on the extracted
Bistatic Radar Cross Section proﬁles and the correlation between them for different types of air-
craft. In order to properly estimate the proﬁles, the ADS-B aircraft location information is adjusted
based on extracted Doppler frequency and then used for Bistatic Radar Cross Section estimation.
The classiﬁcation is based on seven types of aircraft grouped by their size into three classes.

Keywords: passive bistatic radar, very high frequency, bistatic Doppler, bistatic radar cross sec-

tion, instantaneous frequency estimation

Preface

The work presented in this dissertation was realized in the Laboratory of Applied Mathematics of
the Department of Mathematics and Physics of Lappeenranta University of Technology, Finland,
between years 2008–2015. During the ﬁrst years of this period the theoretical part of the meth-
ods presented was studied and developed. In July 2012 the developed techniques were tested and
enhanced with use of acquired data on radio signal data. I would like to acknowledge all the insti-
tutions that have provided a ﬁnancial support for carrying out this research, that is, the Department
of Mathematics and Physics of Lappeenranta University of Technology, the Research Foundation of
Lappeenranta University of Technology. This work would not be possible without a information on
radio signal data, as well as the ADS-B Mode S protocol message data. I would like to thank to my
dear colleagues Juha Hartikka and Mauno Ritola for recording and sharing with me data on radio
signal. Moreover I would like to acknowledge representatives of ﬂightradar24.com for providing
ADS-B data free of charge. I also acknowledge the Finnish Defence Forces Logistics Command for
their insightful analysis, comments and signiﬁcant suggestions that helped building a more reliable
system.

As the research progressed many people had inﬂuence on the direction and momentum of my work.
The ﬁrst and foremost inﬂuential person was my supervisor Tuomo Kauranne. His scientiﬁc guid-
ance has been a cornerstone of this work. Very often the suggestions given by him described the
big picture rather than the details, which helped me with pursuing my goals, inspiring me to think
outside-the-box and most signiﬁcantly straightening some curvy and challenging research paths that
I have encountered. Moreover the abstract way of sharing his thoughts has been clear and well re-
ceived, which is highly desirable in mathematicians’ communities. Tuomo Kauranne has been a
mentor to me for the last eleven years, helping me with both professional and private live, facing
problems together like with a member of the family.

I would like to acknowledge the reviewers of this work, Pekka Neittaanmäki and Hannu-Heikki
Puupponen for their comments which helped creating the ﬁnal form of this thesis. The comments I
have received were gratefully appreciated, being straightforward and to the point.

I would also like to thank Matti Heiliö for his help with realizing the projects related to education,
paper industry, web development. Most of all I want to thank him for giving me a position of
assistant editor of ECMI Newsletter for years 2007–2014 as well as for guiding me through the ﬁrst
years of my stay in Finland. Another person who has helped me understanding the scientiﬁc world
was Heikki Haario. The ﬁrst project that I have realized was supervised by him and was related to
Fast Fourier Transform application for Accurate Period Estimation of Periodic Signals.

I have also received many suggestions outside of an academy, such as DX-er community sharing
their knowledge on electromagnetic signal propagation, aircraft scattering, data reception techniques
and many other topics. I greatly appreciate their help on this matter.

My eleven years lasting stay in Finland was accompanied substantially by Finnish friendship. First,
I want to thank to Miika Tolonen and Anne Tolonen for their hospitality, many times taking me over
to their place and having a great time, showing me what the real Finnish sauna is and how important
cultural role it plays in life, undoubtable mutual trust and ﬁnally for being my friends. It is indeed
a great experience and pure feeling in its nature to become a real friend to Finn, and I am proud to

ONTENTS

Abstract

Preface

Contents

List of the original articles and the author’s contribution

Part I: Overview of the thesis

Introduction

1.1

Review of aircraft tracking systems based on on-board Global Positioning System .

1.2

Motivation for the thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1.3

Non-cooperative system

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

1.4

Radars . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1.5

Bistatic radar . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1.5.1

Passive bistatic radar . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1.5.2

History in brief . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1.5.3

Nonmilitary applications . . . . . . . . . . . . . . . . . . . . . . . . . . .

Used techniques

2.1

Discrete Fourier Transform . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2.2

Short Time Fourier Transform . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2.3

Cell Averaging Constant False Alarm Rate . . . . . . . . . . . . . . . . . . . . . .

2.4

Vincenty’s Inverse Formulae . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2.5

Canny Operator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2.6

Hough Transform . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2.7

Bistatic radar equation

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

2.8

Bistatic Radar Cross Section . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Experiments

3.1

Passive bistatic radar scenario

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

3.2

Transmitter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3.3

Receivers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3.3.1

Dipole . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3.3.2

Yagi-Uda . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3.4

Radio signal data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3.5

Mode S data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Long-distance passive multi-static aircraft tracking

4.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4.2

Previous research on bi and multi-static passive radar . . . . . . . . . . . . . . . .

4.3

Mathematical model

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

4.3.1

Preprocessing the data . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4.3.2

The model

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

4.4

Processing of an Experimental Data Set . . . . . . . . . . . . . . . . . . . . . . .

4.5

Case study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4.6

Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Instantaneous Doppler signature extraction

5.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5.2

PDF of FODDS with respect to varying sampling time and cruising velocity . . . .

5.3

A Doppler curve detection model based on PDF of FODDS . . . . . . . . . . . . .

5.3.1

Cell Averaging – Constant false alarm rate

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

5.3.2

Grouping . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5.3.3

Center of mass . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5.3.4

Expected value . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5.3.5

Classiﬁcation and prediction . . . . . . . . . . . . . . . . . . . . . . . . .

5.3.6

Intersection of sequences . . . . . . . . . . . . . . . . . . . . . . . . . . .

5.3.7

Combining sequences

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

5.4

Data set speciﬁcation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5.4.1

Recorded sessions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5.4.2

Simulated signal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5.5

Case Study

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

5.5.1

Recorded sessions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5.5.2

Simulated signal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5.6

Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Aircraft classiﬁcation based on bistatic radar cross section

103

6.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103

6.2

Data acquisition and preprocessing . . . . . . . . . . . . . . . . . . . . . . . . . . 104
6.2.1

Acquisition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104

6.2.2

Preprocessing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105

6.3

Bistatic radar cross section comparison . . . . . . . . . . . . . . . . . . . . . . . . 109

6.4

Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112

Bibliography

113

OMENCLATURE

Effective area of antenna. 37

Number of properly extracted Doppler signatures. 82, 83, 85, 86

Number of visible Doppler signatures. 82, 83, 85, 86

E/N

Received energy to receiver noise spectral density required for detection. 36

(n) Efﬁciency in adding n

pulses. 37

Preamp noise of the receiver. 81

F Discrete Fourier Transform operator. 24

Accumulation of propagation effects. 37

Receiver noise ﬁgure. 36

Pattern propagation factor for aircraft to receiver path. 36

Pattern propagation factor for transmitter to aircraft path. 36

G Hop size (in samples) between successive DFTs (STFT). 26, 48, 53, 79, 92, 95

Receiving antenna power gain. 36, 81, 97

Transmitting antenna power gain. 36, 81, 97

Gaussian function. 32, 34

Transmitting/receiving antenna gain. 37

The n-order modiﬁed Bessel function of the third kind. 39

The n-order modiﬁed Bessel function of the ﬁrst kind. 27, 39

L Length of window function (STFT). 26, 28, 48, 53, 79, 92

Receiver system losses (>1). 36

Free space path loss of the signal between the transmitter and the receiver. 97

Transmitter system losses (>1). 36

Glossary

Difference in longitude (VIF). 32

Length of the dipole. 46

Length of ogive. 81

System losses. 37

N Number of consecutive integer (positive) numbers. 24

Associated Legendre function of order n and degree 1. 39

Transmitted average power. 36, 80, 85, 86, 97

RC Reference cells (CA-CFAR). 28, 74

Mean radius of the Earth. 94

T Sampling time. 24, 36, 69–72, 75, 76, 82

Standard temperature 290 K. 36

Recording period. 53

Session duration. 83

T r Canny operator threshold. 33

1,2

Reduced latitude (VIF). 32

U b

∗

Set of points on transmitter - receivers baselines. 56

Average cruise speed. 54, 55, 69–72, 75, 76, 80, 85

c,max

Maximum cruising velocity. 96

W Smoothed image (Canny edge operator). 32, 33

The n-order modiﬁed Bessel function of the second kind. 39

α Azimuth of the geodesic at the equator (VIF). 32

Azimuths of the geodesic (VIF). 32

β Bistatic angle. 19, 20, 37, 81, 86, 92, 97

◦ Hadamard product. 76

δ Aspect angle. 19, 20, 37, 92

Pulse width. 81

η Parameter of Kaiser window function. 27

Glossary

γ Angle between the receiver-transmitter vector and the vector of the aircraft’s trajectory measured

counterclockwise. 69, 70, 72

λ Wavelength. 36, 45, 81, 97

(S/N )

Signal to noise ratio with one pulse present. 37

Z Set of integer numbers. 77

Standarized and simpliﬁed form of matrix EC

l1,l2

. 76

l1,l2

Matrix of energy concentration parameters. 75, 76

Simpliﬁed form of matrix F

m
l1,l2

. 76

Standarized and simpliﬁed form of matrix F

l1,l2

. 76

l1,l2

Matrix of parameters that checks every pair of newly found pretenders and pretenders from

the previous scan for their frequency differences. 75, 76

m
l1,l2

Matrix of constrained F

l1,l2

and its boolean values. 75

H Group (set) of points. 74, 75

M Measure of the quality of matching between groups from the previous scan and those from the

present one (matrix). 73, 76, 77

P Simpliﬁed form of matrix P

l1,l2

. 76

l1,l2

Matrix based on PDF of FODDS. 76

S Spectrogram matrix. 26, 28, 32, 53–55, 73–75, 84

F Functional operator that converts arbitrary function to its Discrete Fourier Transform. 24

A Aircraft. 19, 20, 80, 85, 94

E Expected value. 72, 75

FA Number of false alarms. 82, 83, 85, 86

Receiver J

. 43, 44, 48, 78–80, 84, 92, 95

Receiver J

. 43, 44, 48, 58, 78–80, 84

J Receiver J. 43, 44, 50, 53, 55, 58, 59, 61, 63–66, 69, 79, 80, 93

M Receiver M. 43, 44, 48, 50, 53, 55, 58, 59, 61, 63–65, 69, 78–80

Pr Probability density function. 70–72, 76

R Receiver. 18, 19, 21, 31, 32, 43, 70, 92–94

T Transmitter. 18–21, 31, 32, 43, 44, 50, 53, 55, 57–59, 63, 64, 70, 79, 80, 92–94

Glossary

ω Frequency domain. 28, 30, 73–75, 80, 88, 95

Extracted carrier frequency. 93, 95, 96

Frequency of the center for a group. 74, 78, 85, 88, 93, 95, 96

Frequency of the center for a group. 78

Average time span between Doppler signatures. 82, 83, 85, 86

ψ Difference in longitude on an auxiliary sphere (VIF). 32

ρ Distance from the origin to the line along vector perpendicular to the line (HT). 34, 35, 54

Correlation parameter between twosignals. 98

σ Standard deviation. 32

Standard deviation of differences between the smoothed line f

and the Hough Transform (HT)

family of lines. 55

Bistatic radar cross section parameter. 36, 37, 81, 97

Bistatic radar cross section parameter for case of forward scatter. 38

Monostatic radar cross section. 37

Bistatic radar cross section for electric ﬁeld. 38

Standard deviation of difference f

− ω

. 85–88

τ Discrete time. 24–26, 53

θ The angle the normal line makes with x-axis (HT). 34, 35, 54, 60

ϕ Half angle of ogive’s nose. 81, 86

ξ angular distance on the auxiliary sphere from T to R (VIF). 32

angular distance on the sphere from the equator to T (VIF). 32

angular distance on the sphere from the equator to the midpoint of the line (VIF). 32

Great circle arcs connecting the aircraft A with the receiver R. 94

Great circle arcs connecting the aircraft A with the receivers J and M. 57

Great circle arcs connecting the transmitter T with the aircraft A. 57, 94

a Major semiaxis of the ellipsoid (VIF). 32

aSNR Average signal to noise ratio. 82, 83, 85–88

Extracted amplitude of the carrier. 93, 96–98, 100

Glossary

Amplitude of the center for group H. 74, 85, 93, 95–98, 100

Amplitude of synthetic Doppler signal. 85

alt Altitude of the aircraft A. 57, 80, 85, 94, 96

b Minor semiaxis of the ellipsoid (VIF). 32

Cardinality of group (set) H. 74, 75

Baseline for receiver J. 56, 58

Baseline for receiver M. 56, 58

c Velocity of propagation of electromagnetic waves (light). 56, 69, 94, 96

Distance between receiver J and receiver M. 43, 44, 79, 80

T(R)A

Monostatic transmitter(receiver) to the aircraft distance. 36, 37

Distance between transmitter T and receiver J. 43, 44, 79, 80

Distance between transmitter T and receiver M. 43, 44, 79, 80

Bistatic distance. 69

Distance from aircraft A to receiver R. 36, 56, 94, 97

Distance from aircraft A to receiver I = J, M. 69

Distance from transmitter T to aircarft A. 36, 56, 69, 94, 97

Baseline length. 19, 31, 32, 69, 70, 80, 92, 94, 97

Baseline length. 97

min

Maximal minimum distance that have to be attained between the trajectories. 98

ec Energy concentration value. 75

en Scan numbers when the sequence ended. 76–78, 85, 87, 88

f Frequency. 24

Smoothed curve. 55

Discontiuous curve based on f

and HT family of lines ﬁt. 55

Doppler frequency. 56, 69–72, 76, 85, 88

FR24

Second estimation of Doppler frequency based on FR24 location/time data. 95, 96

FR24

Doppler frequency based on FR24 location/time data. 94, 95

Flattering of the ellipsoid (VIF). 32

Glossary

Horizontal ﬁlter. 32

Vertical ﬁlter. 32

Doppler curve that corresponds to the best ﬁt. 58, 62

Curve based on f

with outliers removed and replaced with interpolated values. 55–57, 62

Frequency margin. 54, 59

Pulse repetition frequency parameter. 37, 81

Sampling frequency. 25, 28, 53, 81

Transmitting frequency for transmitter T. 43, 46, 48, 53, 54, 59, 69, 79–81, 92, 94, 96

D,max

Maximum achievable Doppler shift. 95, 96

ar1

Discrete signal (function dependent on discrete time). 25

ar2

Discrete signal (function dependent on discrete time). 25

Discretized form of h

. 24, 26

cost

Cost function. 95

Frequency deviation limit for sequence w

to be categorized as a carrier. 77, 82

mrp

Frequency margin for predicted values. 78

Frequency margin for grouping procedure. 74, 75, 82

Quality measure for a sequence. 76, 77

Arbitrary contiuous function. 24

j Frequency index. 28, 73–75

k Boltzmann’s constant. 36

Wave number. 39

CF AR

Constant for CA-CFAR. 28, 74, 82

Anticipated points. 62, 64

lat Latitude. 94, 97

lat Estimated latitude. 96, 100

lat

Shift in latitude. 95, 96

lon Longitude. 94

lon Estimated longitude. 96, 100

Glossary

lon

Shift in longitude. 95, 96

m Size of the spectrogram matrix in frequency domain. 28, 73, 74

n Size of the spectrogram matrix in time domain. 28, 53, 73

Length of the sequence for which we check if its trend (ﬁrst order polynomial ﬁt) p

has deviated

by f

. 77, 82

Length of the signal that is needed for predicting its next frequency value. 75, 76, 82

Length of sequence for which sequence is considered as a signal. 76–78, 82

Size of matrix S in the time interval between t

and t

. 54

Length of reference cells (CA-CFAR). 28, 74, 82

Upper limit for parameter n

. 76, 82

l,z

Length of prediction to the past and to the future for w

and w

respectively. 78

Number of pulses. 37

p Time index. 28, 30, 69, 73–77

First order polynomial ﬁtted to w

. 77, 82

Termination coefﬁcient. 77, 82

q Number of pairs (ω

) (pretenders). 74–76

Ratio between lengths of extraction time t

and exact Doppler time t. 85–88

Radius of the sphere. 39

st Scan numbers when the sequence started. 76–78, 85, 87, 88

t Continuous time. 24, 28, 30, 69–78, 85, 87, 88, 93–96

Time occurence of the ﬁrst sample. 24, 53, 63–65

Calculation time needed for tracing the spectrogram image. 83, 85, 86

Lower time limit where the Doppler curve was found. 54, 60, 61

Time margin for the potential signal. 78, 82

Upper time limit where the Doppler curve was found. 54, 60, 61

Crossing time. 55, 60, 61

Second estimation of crossing time. 55, 56, 61

Aircraft’s trajectrory i. 97, 98

1.3 Non-cooperative system

1.3

Non-cooperative system

A Non-Cooperative System (NCS) is a system that does not require any information from the aircraft
about its state, even if such information could be provided by the aircraft. In this case we consider
the aircraft as a Non-Cooperative Target (NCT).

An example of such a system is presented in this work as Passive Bistatic Radar (PBR) conﬁguration
in which case the means of tracking that are independent on any collaboration by the aircraft are
used.

1.4

Radars

The deﬁnition of radar as formulated in IEEE (1990) states:

Deﬁnition 1.4.1 A device for transmitting electromagnetic signals and receiving echoes from ob-
jects of interest (targets) within its volume of coverage. Presence of a target is revealed by detection
of its echo or its transponder reply. Additional information about a target provided by a radar in-
cludes one or more of the following: distance (range), by the elapsed time between transmission of
the signal and reception of the return signal; direction, by use of directive antenna patterns; rate of
change of range, by measurement of Doppler shift; description or classiﬁcation of target, by anal-
ysis of echoes and their variation with time. The term radar was originally an acronym for radio
detection and ranging.

According to Farina (2005) radar systems can be categorized by its features into the following
subsets:

i. Radar location: ground-based: ﬁxed, transportable, mobile; ship-borne; air-borne; space-

borne;

ii. Capacity: tracking, surveillance, reconnaissance, imaging;

iii. Applicability:

• air defence,

• air trafﬁc control (terminal area, en route, collision avoidance, apron),

• monitoring of surface trafﬁc in the airports (taxi radar),

• anti ballistic missile defence,

• vessel trafﬁc surveillance,

• remote sensing (application to crop evaluation, hydrology, geodesy, archaeology, astron-

omy, defence),

• meteorology (hydrology, rain/hail measurement),

• study of atmosphere (detection of micro-burst and gust, wind proﬁlers),

• space-borne altimetry for measurement of sea surface height,

• acquisition and tracking of satellites in the re-entry phase,

1. Introduction

• monitoring of space debris,

• anti-collision for cars,

• ground penetrating radar (geology, gas pipe detection, archaeology, detection and loca-

tion of mines, etc.);

iv. Band (see IEEE Standard (521

) (2003));

v. Beam scanning:

• ﬁxed beam,

• mechanical scan (rotating, oscillating),

• mechanical scan in azimuth,

• electronic scan (phase control, frequency control and mixed in azimuth/elevation),

• mixed (electronic-mechanical) scan,

• multi-beam conﬁguration;

vi. Number and type of collected data:

• range (delay time of echo),

• azimuth (beam pointing of antenna beam, amplitude of echoes),

• elevation (only for three-dimensional radar, multifunctional, tracking),

• height (derived by range and elevation),

• intensity (echo power),

• Radar Cross Section (RCS) (derived by echo intensity and range),

• radial speed (measurement of differential phase along the time on target due to the

Doppler effect; it requires a coherent radar),

• polarimetry (phase and amplitude of echo in the polarisation channels: HH - horizontally

transmitted, horizontally received - HV, VH, VV),

• RCS proﬁles along range and azimuth (high resolution along range, imaging radar);

vii. Conﬁguration:

• monostatic (co-located T and R - same antenna, mono-radar/multi-radar),

• bistatic (not co-located T and R - two antennas),

• multistatic (one or more T and R spatially dispersed);

• suitable references for bistatic, multistatic and passive radar are: [2], [19] to [21];

viii. Waveform: Continuous Wave (CW), pulsed wave, digital synthesis;

ix. Processing:

• coherent (MTI/MTD/Pulse-Doppler/super-resolution/SAR/ISAR),

• non coherent (integration of envelope signals, moving window, adaptive threshold (CFAR))

1.5 Bistatic radar

Radar

Bistatic and

multistatic

radar

Coop-

erative

trans-

mitter

Non-

cooperative

trans-

mitter

Radar

transmitter

Broadcast/

comms/

radion-

avigation

transmitter

Monostatic

and quasi-

monostatic

radar

Figure 1.1: Bistatic and multistatic radar’s taxonomy.

• mixed;

x. Technologies:

• for antenna (reﬂector plus feed, array (planar, conformal), corporate feed/air - cou-

pled/lens),

• transmitter (magnetron, klystron, TWT, mini TWT, solid state)

• receiver (analogue and digital technologies, base band, intermediate frequency sam-

pling, etc.; relevant parameters of receiver are: noise ﬁgure, bandwidth and dynamic
range).

The bold text in the above list points to the characteristics of the radar that were used in the further
chapters.

Additionally Grifﬁths (2010) attempts to categorize the radar family into the subsets as shown in
Fig. 1.1.

1.5

Bistatic radar

Deﬁnition of BR is stated in (IEEE, 1990):

Deﬁnition 1.5.1 A radar using antennas at different locations for transmission and reception.

In Fig. 1.2 transmitter T and receiver R are separated by a line of length L, in contrast to monostatic
radar where both sides are colocated, which is called a baseline. In further studies we will denote
the baseline as d

. The target is denoted by A and in this case it is an aircraft. In general any

1. Introduction

of these items may be categorized as ground-based, airborne or marine and can be moving or be
stationary. The angle between vectors deﬁned by the illumination path and the echo path (positive
angle, less than 180

◦

) is called a bistatic angle and is denoted with β. It is also sometimes called

the cut angle or the scattering angle. Another crucial angle is the aspect angle denoted by δ and it
deﬁnes the target’s movement at speed v with respect to the bistatic bisector.

Direct path

Baseline L = d

Echo

path

Illumination

path

β/2

Figure 1.2: Bistatic triangle with accompanied features and variables.

The distance between the transmitter T and the receiver R is not explicitly given in the deﬁnition
of bistatic radar. However, in Skolnik (2003) the author quantiﬁes the separation as one of “a
considerable distance” or “comparable with the target distance” or “the echo signal does not travel
over the same path as the transmitted signal”.

1.5.1

Passive bistatic radar

One of the forms of Bistatic Radar is the Passive Bistatic Radar. Passive Bistatic Radar is also
known as bistatic hitchhiker or parasitic radar and it uses illuminators of opportunity as transmitters
to track an object.

In general the concept of the PBR can be speciﬁed as:

• PBR is a subtype of BR (all bistatic analysis such as geometry, doppler, RCS apply),

• PBR is a BR that does not emit any Radio Frequency (RF) of its own to track targets,

• it utilizes the existing RF energy in the atmosphere,

• as sources of RF energy we can use broadcast Frequency Modulation (FM) stations, GPS,

cellular telephones and commercial television,

1.5 Bistatic radar

Merits

Disadvantages

no demand for frequency allocation

no direct control over emitted signal

relatively low cost

complicated geometry with respect to
monostatic radar

covert receiver location, no possibility of
jamming

technology is outdated

immune

Anti-Radiation

Missiles

(ARM)

possibility to detect stealth objects

good low level coverage

large number of transmitters

complements existing system

Table 1.1: Merits and disadvantages of PBR usage.

• Hitchhiker or parasitic radar are used in a case when another radar transmitter is used, such

as signal from monostatic radar,

• when the transmitter of opportunity is from a non-radar transmission, such as broadcast com-

munications, then terms such as passive radar, passive coherent location or passive bistatic
radar are used.

The common use of PBR is to use RF energy such as commercial FM broadcast as a transmitter T
which is scattered by an aircraft A. Reception of the scattered signal is conducted with an antenna
and compared with the reference signal (direct path signal) from a second receiving antenna. Then
by using Digital Signal Processing (DSP) techniques, target-related parameters such as range, range-
rate, Angle of Arrival (AoA) and other can be estimated.

An alternative method is the Doppler-only technique which does not use a second receiving antenna,
but rather a system of bistatic radars working in unison. The received signal is analyzed by putting
great importance on Doppler shift information. The information acquired by the receiving parties
involved is then combined in order to detect/track targets.

Merits and disadvantages of using a PBR conﬁguration are presented in Table 1.1.

Examples of such transmitters are analog radio/tv transmitters, cellular phone base stations and
Digital Audio Broadcast (DAB) that have been presented in Grifﬁths and Baker (2005); Malanowski
et al. (2014) or Global Navigation Satellite System (GNSS) in Clemente and Soraghan (2014);
Suberviola et al. (2012).

1.5.2

History in brief

First experiments on Bistatic Radar have been conducted simultaneously in United Kingdom, The
United States, France, the Soviet Union, Japan, Germany and Italy in the late 1930s, before the
Second World War. Some of the experiments were deployed as forward scatter fences (along country

1. Introduction

Figure 1.3: A Klein Heidelberg in Tausendfüssler near Cherbourg (left) and at Biber (middle).
Towers of British Chain Home (right).

borders) or as a bistatic hitchhiker as an aircraft detection systems. These radars were known as
Continuous Wave detectors and served to detect an object as it crossed the baseline by estimating
the frequency of the Doppler shift caused by an airborne object with respect to the direct signal
path from transmitter T to a parasitic or dedicated receiver R. In 1930 an aircraft was accidently
spotted by L.A. Hyland while working on the directional properties of an aircraft antenna system at
32.8 MHz and being located 3.2 km from the aircraft. Two years later a team of Taylor, Young and
Hyland used bistatic radar to detect aircraft 80 km from the transmitter. In 1935 Britons Watson-
Watt and Arnold Wilkins conducted an experiment of aircraft detection in a forward-scatter fence
conﬁguration which later on was developed into the Chain Home network (see Fig. 1.3) of radars
along the British coast. The network appeared to be effective during the Battle of Britain in 1940.

In Germany bistatic radar was developed during World War II under the name Klein Heidelberg, see
Fig. 1.3. It was using the British Chain Home radar as an illuminator of opportunity and was used
to detect approaching allied bombers when they were passing the English Channel.

The ﬁrst wave of interest in BR ended in mid 1930s. The interest was revived in 1950s when the
development in monostatic radar theory, techniques and technology could be used also in a bistatic
conﬁguration. During this time a couple of new theories were introduced including Bistatic Radar
Cross Section (BRCS) and bistatic clutter measurement. Development of semiactive homing mis-
siles and a second attempt at hitchhiking and forward-scatter fences took place. This time hitchhik-
ing was aimed at locating objects in space using atmospheric phenomena, such as lighting or radio
stars like the Sun, as transmitters. Multistatic radars were designed with the purpose of ballistic
missile detection in systems like the U.S. Plato, satellite detection and tracking fence Navsparus.
Only the last one was actually deployed.

The second resurgence of the BR subject is dated to 1970s. During these years several new systems
were introduced including Sanctuary, presented in Bailey et al. (1977), Bistatic Alerting and Cueing
(BAC) in Thompson (1989), Aircraft Security Radar (ASR) and Multistatic Measurement System
(MMS) in Willis (2005).

The main objective of the Sanctuary project was to generate a clutter-free display for air defense
operations with a transmitter onboard of an aircraft and ground-based coherent receiver. Success-
ful ﬂight tests were performed during which jet attack aircraft was detected at range greater than
100 km.

The ASR on the other hand had a purpose of detecting intruders approaching an aircraft. The

2. Used techniques

In this work STFT is used to represent a signal in the time-frequency domain.

When calculating the Fourier Transform (FT) of some given function, it might happen that it is
deﬁned in terms of a continuous independent variable, like in the case of transform pairs. In the
case when function values are given as discrete values at regular time intervals, like with physical
measurements, the transformed function will also be available at discrete intervals. One may think
of a discrete function as a sort of approximation of an underlying function of a continuous variable.

To understand the Discrete Fourier Transform (DFT) one must understand the idea of periodicity.
A periodic function is a function the values of which are repeated at equal intervals of T seconds.
We can say then that the values are repeated with frequency

Hz.

Let us assume that the parameter τ represents a time with N consecutive integer (positive) values.
To illustrate the discretizing operation, let us take the sine function on interval [0,

3
2

π], see Fig. 2.1.

−1

)

−1

[τ

]

Figure 2.1: An example of the discretization of the sine function. Note the change in scale on
horizontal axes for continuous time t and discrete time τ .

The domain notation has changed from continuous time notation t to a discrete one τ , thus for the
function sine, values of h

are only known at discrete time instances f

[τ ]. We can summarize this

operation by formulae 2.1.

[τ ] = h

+ τ T )

(2.1)

where T stands for the sampling interval and t

denotes time occurrence of the ﬁrst sample. The

deﬁnition of DFT given in Bracewell (2000) states that f

[τ ] possesses a discrete Fourier transform

F (v) expressed by

F (v) = N

−1

N −1

τ =0

f [τ ] e

−i2π(v/N )τ

(2.2)

It is important to stress that there are different frequency notations for the continuous case and the
discrete one. The frequency notation for FT is f , whereas for DFT it is v/N which can be interpreted
as quantity measured in cycles per sampling interval.

The DFT shares all the properties of FT. However the form of these properties may be different.
Using the functional operator F , which converts a function to its Fourier transform, we assume that

2.1 Discrete Fourier Transform

F [f

ar1

[τ ]] = F

ar1

iτ ω

) and F [f

ar2

[τ ]] = F

ar2

iτ ω

), where ω = 2πv/N , f

ar1

and f

ar2

are two

discrete signals (functions dependent on discrete time τ ).

• linearity

F [af

ar1

[τ ] + bf

ar2

[τ ]] = aF

ar1

iτ ω

+ bF

ar2

iτ ω

• time shifting

F [f

ar1

[τ − τ

]] = e

−iτ

ar1

iτ ω

• time reversal

F [f

ar1

[−τ ]] = F

ar1

−iτ ω

• frequency shifting

F f

ar1

[τ ] e

iω

= F

ar1

iτ (ω−ω

)

• differencing

F [f

ar1

[τ ] − f

ar1

[τ − 1]] = 1 − e

−iτ ω

ar1

iτ ω

• differentiation in frequency

−1

dω

ar1

iτ ω

= τ f

ar1

[τ ] ,

• convolution theorems

F [f

ar1

[τ ] ∗ f

ar2

[τ ]] = F

ar1

iτ ω

ar2

iτ ω

F [f

ar1

[τ ] f

ar2

[τ ]] = F

ar1

iτ ω

∗ F

ar2

iτ ω

• Parseval’s relation

∞

τ =−∞

ar1

[τ ]|

2π

ar1

iτ ω

dω

An example of DFT in a form of amplitude signal is presented in Fig. 2.2. Here as an input
signal the summation of two sine functions is taken, such as f

[τ ] = 5 sin (2πτ f

+ π/4) +

2 sin (2πτ f

− π/2) + e, where f

= 45 Hz, f

= 110 Hz and e ∼ N (0, 1) denotes a source of noise

that follows the normal distribution. Sampling frequency is equal f

= 1 kHz which means that time

length of the sample presented in Fig. 2.2 is equal τ /f

= {1, ..., 500} /1000 = {0.001, ..., 0.5} s,

0.5 s.

The spectrum in Fig. 2.2 depicts the signal after the transformation. The peaks represent the main
periodic components of the signal under consideration. The ﬁrst peak is located at 44.92 Hz, the sec-
ond one at 109.4 Hz. This inaccuracy in frequency determination compared to the original settings
is due to the length of the signal – the more samples of the signal the better the frequency estimation.
The reason the amplitudes on the spectrum are not exactly 5 and 2 (these are estimated as 5.04 and
1.74, respectively) is due to the noise e. The remedy for this is as in the case of frequency inaccuracy
is to increase the length of the considered signal.

2. Used techniques

100

200

300

400

500

−5

discrete time τ

[τ

]

100

200

300

400

500

frequency Hz

Figure 2.2: Original discretized signal (left) and its single-sided amplitude spectrum (right)

2.2

Short Time Fourier Transform

The STFT is a signiﬁcant tool of time-frequency representation in ﬁelds like radar (Rothwell et al.,
1998; Pan et al., 2013), sonar (Ferguson, 1996), music (Pielemeier et al., 1996), speech (Liu, 1993),
(instantaneous) Frequency Modulation (FM) demodulation (Kwok and Jones, 2000). Due to the
demand of quick execution time, an efﬁcient instantaneous implementation of the STFT is essential.
The merits of STFT over other time-frequency distributions are presented in (Durak and Arikan,
2003).

In order to transform a signal f

with STFT (Allen and Rabiner, 1977), also known as short-time

spectrum (Hlawatsch and Boudreaux-Bartels, 1992), into matrix form (spectrogram image) with
time-frequency (t, ω), domain expressed in s and Hz, the following equation 2.3 must be used

S (t, ω) =

N −1

τ =0

[τ ] w [τ − tG] e

−iωτ

(2.3)

where f

(τ ) is an input signal at time τ , w is a length L window function, S(t, ω) is the Discrete

Fourier Transform of windowed data centered about time tG and G denotes a hop size (in samples)
between successive DFTs. Further we will refer to S, S(t, ω) as the spectrogram resulting from the
STFT operation.

The procedure of applying STFT to the signal might be outlined as follows:

• Split the signal into blocks of equal length L,

• The blocks overlap each other by some factor 100 % > 1 −

% ≥ 0 %,

• Each consecutive block is windowed. The windowing operation is a multiplication of the

block with a smooth function with tails that are nearly zero.

In Oppenheim and Schafer (2009), the authors list window functions typically used for smoothing:

2.2 Short Time Fourier Transform

• Rectangular

w [n] = 1, 0 ≤ n ≤ L,

• Hann

w [n] = 0.5 − 0.5 cos

2πn

, 0 ≤ n ≤ L,

• Hamming

w [n] = 0.54 − 0.46 cos

2πn

, 0 ≤ n ≤ L,

• Blackman

w [n] = 0.42 − 0.5 cos

2πn

+ 0.08 cos

4πn

, 0 ≤ n ≤ L,

• Kaiser

w [n] =

η 1 −

− 1

1
2

(η)

, 0 ≤ n ≤ L.

with J

, n = 0 the zeroth-order modiﬁed Bessel function of the ﬁrst kind. Rectangular for

η = 0, gaussian for η → ∞.

Shapes of the listed window functions are presented in Fig. 2.3.

0.5

]

Rectangular

Hann

Hamming

Blackman

0.5

]

η = 1

η = 2

η = 5

η = 10

η = 20

η = 50

Figure 2.3: Window functions. Rectangular, Hann, Hamming, Blackman (left). Kaiser with param-
eter η = {1, 2, 5, 10, 20, 50} (right).

One of the signiﬁcant properties of STFT is its ability to represent time-frequency content of signals
that is free of cross terms as presented in Durak and Arikan (2003). The existence of cross-terms
however can be observed when using other transforms such as Wigner Distribution (WD). Here the
cross-term existence is due to the auto-correlation function inherent in its formulation. Cross-term
elimination is an important quality for distinguishing real components from artifacts which is crucial
for the case presented in this work, therefore further analyzes here are based on STFT application.

The other important feature of STFT is its computational simplicity. In general the properties of
STFT with respect to other time-frequency analysis transforms can be summarized as

2. Used techniques

• the clarity of representation is of low quality,

• the ﬁxed resolution issue. The width of a window is a trade-off between time and frequency

resolution. Longer window (narrowband) lead to higher frequency resolution but decreases
time resolution which can be observed as smudges (smoothing) in time direction. The short
window (wideband) determines high time resolution and low resolution in frequency, which
is also undesired,

• cross-term free transform,

• low computational complexity.

By multiplying one of the aforementioned window functions with the window of the signal we
attenuate the middle part and suppress the tails of the windowed signal which characterizes STFT
as a local spectrum of the signal around the analysis window.

An example of the application of Short Time Fourier Transform is depicted in Fig. 2.4.

Here we consider recorded speech of the sentence The quick brown fox jumped over the lazy dog.
The example ought to show how the one dimensional signal (bottom) can be represented in a more
informative and intuitive way which is a two dimensional time-frequency domain (center and top).
The top two ﬁrst images are presented to visualize a difference between different window sizes L.
The ﬁrst image represents STFT with window width L = 1 ms which can be interpreted as a good
trade-off between the resolutions of the signal of 44 100 Hz sampling frequency f

. The purpose

of the second image is to show how the undesirable smoothing factor, lack of sharpness in time
direction, for STFT come into play for window width L = 10 ms.

2.3

Cell Averaging Constant False Alarm Rate

As mentioned in 2.2 a matrix of spectrogram data is denoted by S

[n×m]

and an amplitude of each cell

of this matrix by S (t (p) , ω (j)), where t (p) , p ∈ [1, n] denotes a time for the cell being measured,
and ω (j) , j ∈ [1, m] the frequency that corresponds to the cell. n and m are sizes in time and
frequency domains, respectively, of the spectrogram.

Constant False Alarm Rate (CFAR) technique is well documented, and many variants of it have
been developed.

As an example of this technique we want to give a brief presentation of one of the aforementioned
models, namely the Cell Averaging – Constant False Alarm Rate (CA-CFAR).

Let us consider Fig. 2.5 in which there are three distinguished groups of cells: reference cells
(dotted), guard cells (crosshatch), and the Cell Under Test (CUT) (grid). To check if detection is
declared in the CUT we need to average over all the reference cells RC of length n

and then after

multiplying it with the constant k

CF AR

compare with the amplitude of CUT, see (2.4).

S (t (k) , ω (p)) > k

CF AR

j∈RC

S (t (k) , f (j))

(2.4)

where RC denotes a reference cell’s location in the frequency domain and is of length n

, k

CF AR

stands for a scaling constant. If the inequality is satisﬁed then CUT is stored and denoted as
S (t (k) , ω (p

, t (k))) , p

∈ [1, m].

2.4 Vincenty’s Inverse Formulae

frequenc

kHz

frequenc

kHz

0.5

1.5

2.5

3.5

4.5

time s

The

quick

brown

fox

jumped

over

the

lazy

dog

Figure 2.4: Recording of the sentence The quick brown fox jumped over the lazy dog with sampling
frequency f

= 44 100 Hz (bottom); spectrum representation with Hamming window of length

L = 1 ms and overlapping ratio of 50 % (top); the spectrum with L = 10 ms, 50 % (center).

2.4

Vincenty’s Inverse Formulae

Assuming that ﬂight trajectories are presented in spherical coordinates the following formulae is
applied in order to calculate the geodesics.

In Vincenty (1975), the author presents a new method for deriving the shortest distance between
two arbitrary locations (length of geodesics) on the Earth. The method features a relatively high

2. Used techniques

t (p − 1)

t (p)

t (p + 1)

(1)

)

Reference cells

Guard cells

CUT

Figure 2.5: Cell averaging constant false alarm rate scheme.

precision (down to 0.6 mm accuracy) and low time consumption, so it is used often in applications.
The compactness of the formulae is due to nested equations to compute elliptic terms. From two
techniques presented by Vincenty, the direct and the inverse, we need to use Vincenty’s Inverse
Formulae (VIF), which is an iterative approach and is sketched in the following blocks:

ψ = L

(2.5)

sin

ξ = (cos U

sin ψ)

+ (cos U

sin U

− sin U

cos U

cos ψ)

(2.6)

cos ξ = sin U

sin U

+ cos U

cos U

cos ψ

(2.7)

tan ξ =

sin ξ

cos ξ

(2.8)

sin α = cos U

cos U

sin ψ

sin ξ

(2.9)

cos 2ξ

= cos ξ − 2 sin U

sin U

cos

(2.10)

cos

α (a

− b

)

(2.11)

parameter ψ is derived from equations 2.12 and 2.13

C =

cos

α 4 + f

4 − 3 cos

(2.12)

ψ = L

+ (1 − C) f

sin α ξ + C sin ξ cos 2ξ

+ C cos ξ −1 + 2 cos

2ξ

(2.13)

The block of equations starting from equation 2.6 is being derived until the change in parameter ψ
is negligible (change of about 10

−12

). After that we evaluate the following

= bA (ξ − ∆ξ)

(2.14)

2.4 Vincenty’s Inverse Formulae

R α

Figure 2.6: Sphere with parameters used in Vincenty’s inverse formulae

where ∆ξ can be calculated from equations 2.15, 2.16 and 2.17

A = 1 +

16384

4096 + u

−768 + u

320 − 175u

(2.15)

B =

1024

256 + u

−128 + u

74 − 47u

(2.16)

∆ξ = B sin ξ cos 2ξ

B cos ξ −1 + 2 cos

2ξ

−

B cos 2ξ

−3 + 4 sin

−3 + 4 cos

2ξ

(2.17)

The azimuths may be estimated from equations 2.18 and 2.19

tan α

cos U

sin ψ

cos U

sin U

− sin U

cos U

cos ψ

(2.18)

tan α

cos U

sin ψ

− sin U

cos U

+ cos U

sin U

cos ψ

(2.19)

(2.20)

Notation for the Vincenty formulae.

2. Used techniques

a, b

major and minor semiaxes of the ellipsoid, a = 6 378 137 m, b =
6 356 752.314 245 m

a−b

ﬂattening, f

= 1/298.257223563

difference in longitude

length of the geodesic between T and R

, α

azimuths of the geodesic, clockwise from north (forward az-
imuths). α

is produced in the direction from T to R. α

is pro-

duced in the direction from R to Z

azimuth of the geodesic at the equator

1,2

reduced latitude, deﬁned as tan U

1,2

= (1 − f ) tan ψ

difference in longitude on an auxiliary sphere

angular distance on the auxiliary sphere from T to R

angular distance on the sphere from the equator to T

angular distance on the sphere from the equator to the midpoint
of the line M

Values for a, b and f

are taken from World Geodetic System 84 (WGS84).

2.5

Canny Operator

Canny operator also known as Canny edge detector (Russell and Norvig, 1995) is a standard al-
gorithm used for detecting edges in images, in this case it is used to estimate edges within the
spectrogram S (t, ω). To detect an edge within the spectrogram S at any orientation, the spectro-
gram must be convolved with two ﬁlters f

= G

(x) G

(y) and f

= G

(y) G

(x), where f

and f

are vertical and horizontal ﬁlters, respectively, and G

(x) is a Gaussian function expressed

in equation 2.21

(x) =

√

2πσ

−x2

2σ2

(2.21)

where σ denotes standard deviation. The differentiated form of the Gaussian function is derived in
equation 2.22

(x) =

−x

√

2πσ

−x2

2σ2

(2.22)

The procedure used in Canny edge detector can be speciﬁed as follows

i. Calculate convolution of the spectrogram S with ﬁlters f

(x, y) and f

(x, y). The resulting

images are denoted with W

(x, y) and W

(x, y). Let us also deﬁne W (x, y) = W

(x, y) +

(x, y)

2.6 Hough Transform

ii. Calculate the absolute value of W (x, y)

iii. Search for the values in |W | (x, y) that are higher than some speciﬁed threshold T r.

The marked edge pixels are then linked together to form edge curves. This can be achieved by
making the assumption that any neighboring pixels (cells of matrix W ), that were found after the
threshold operation and that keep the orientation consistent, must belong to the same edge curve.

In Canny (1986) the author summarizes the performance criteria to include the following:

• Performance in detection ensures a low probability level of misdetection of real edge points,

and low probability level of detecting non-edge points. Maximizing Signal to Noise Ratio
(SNR), both probabilities can achieve signiﬁcantly lower values.

• Detected edge points (cells) should be located as close as possible to the real edge’s center.

• Multiple detections of a single point must be avoided. The ﬁrst criterion ensures that this

condition is ﬁlled by rejecting points that are considered false.

A one-dimensional example of applying the Canny operator is presented in Fig. 2.7. The convo-
lution of the signal presented in (i) is shown in c), the ﬁrst derivative of the Gaussian function
presented in equation 2.22 is shown in b). The maximum of the absolute value of the convolved
function represents the edge point, here peaking at around 150.

2.6

Hough Transform

The Hough Transform (HT) is a global processing method for boundary detection. The features that
may be extracted with use of HT are circles, lines, ellipses and two-dimensional shape identiﬁcation
in general. There are many different versions of this technique, some of them aiming at reducing
algorithm’s time consumption and some at decreasing its memory usage like the one presented in
Chutatape and Guo (1999).

Firstly a binary edge image is estimated, with the use of the Canny operator f.e., after which each
edge point is transformed to a parameterized curve. The next step is to accumulate the accumulator
array which is usually implemented as an array of accumulator space. Each image pixel (cell of
matrix W ) gives one score to the cells lying on its transformed curve. The last step estimates the
local maxima. Cells with the local maximum of scores have coordinates of parameter representing
a curve (line) segment in the image space (matrix W ).

The kernel of the standard version of HT can be outlined as follows:

i. form the set W of all edge points in a binary image,

ii. transform each point of the set W into a parameterized curve of the parameter space,

iii. accumulate the accumulator space under each parametric curve,

iv. ﬁnd accumulator cells that contain local maxima corresponding to image space curve param-

eters.

2.7 Bistatic radar equation

the highest which means it represents a high accumulation level. These points are a parametric
representation of the lines. The next step is to identify those local maxima from image space curve
parameters. Then the identiﬁed points are used for searching for line segments within the image
presented in top left. The result of this procedure is depicted in c).

1,000

2,000

3,000

1,000

2,000

a) Original image

1,000

2,000

3,000

1,000

2,000

c) Extracted features

−80

−60

−40

−20

−4,000

−2,000

2,000

4,000

b) HT of the image

Figure 2.8: Hough Transform applied to an image. a) the image with two lines 1 and 2; b) HT of
the image in parameterized space (θ, ρ) with local maxima denoted as 1 and 2 corresponding to the
lines; c) extracted lines.

2.7

Bistatic radar equation

This section is presented assuming that in both cases of monostatic and bistatic radar the radar
waveform (transmitted signal) is in a form of repetitive series of rectangular pulses. Moreover its
purpose is to familiarize the reader with the concept of Bistatic Radar Cross Section (BRCS) which
is deﬁned by the Bistatic Radar Equation (BRE).

Following the deﬁnition described in Skolnik (2008) and in Willis and Grifﬁths (2008) the equation

2. Used techniques

is expressed in (2.24) as

)

max

T G

(4π)

(E/N

) L

1
2

(2.24)

where

transmitter to aircraft distance

[m]

aircraft to receiver distance

[m]

transmitted average power

[J s

−1

]

signal observation time

[s]

transmitting antenna power gain

[1]

receiving antenna power gain

[1]

wavelength

[m]

BRCS parameter

]

pattern propagation factor for transmitter to aircraft path

[1]

pattern propagation factor for aircraft to receiver path

[1]

Boltzmann’s constant

1.38 × 10

−23

[J K

−1

]

standard temperature

290

[K]

receiver noise ﬁgure

[1]

E/N

received energy to receiver noise spectral density required for
detection

[1]

transmitter system losses (>1)

[1]

receiver system losses (>1)

[1]

Comparison of the bistatic range equation with that of monostatic conﬁguration (2.24) shows that the
form of the equation is similar. The deﬁnite difference is the replacement of d

2
T(R)A

with d

where d

T(R)A

denotes monostatic transmitter(receiver) to the aircraft distance. This replacement

causes fundamental changes to the operational system of Bistatic Radar (BR) compared to the
monostatic. One of these changes is that the region where the bistatic radar can operate is now
deﬁned by ovals of Cassini, rather than a circle – as in the case of monostatic radar.

Non-collinearity of bistatic constant range sum, deﬁned by d

+ d

, with bistatic constant de-

tection contours (ovals of Cassini) is another signiﬁcant difference between monostatic and bistatic
conﬁgurations – in the case of monostatic radar these two contours, which are the circles, over-
lap each other. This causes aircraft’s E/N

parameter to change with respect to its location on a

constant range sum contour - in contrast to the constant value of E/N

for monostatic radar case.

The next difference is a varying width of bistatic range cell with respect to the target’s location on a
constant range sum contour (ellipse) – against a constant width for monostatic.

The fourth difference is given in a form of pattern propagation factors F

and F

which in the

case of monostatic radar are usually equal, whereas for the bistatic radar they almost always differ.

These four differences are given so as to understand the operational phenomena of bistatic radar
based on the easy-to-grasp monostatic radar operation. Therefore, for reference purposes, the mono-
static equation is presented in (2.25). Note that in the case of monostatic radar, usually both trans-

2.8 Bistatic Radar Cross Section

mission and reception are managed with the same antenna antenna.

2
T(R)A

(n) F

(4π)

(S/N )

1
2

(2.25)

where

T(R)A

maximum radar range

[m]

antenna gain

[1]

effective area of antenna

]

monostatic radar cross section of the aircraft

]

number of pulses

[1]

(n)

efﬁciency in adding n

pulses

[1]

accumulation of propagation effects

[1]

Pulse Repetition Frequency (PRF) parameter

[Hz]

(S/N )

signal to noise ratio with one pulse present

[1]

system losses

[1]

2.8

Bistatic Radar Cross Section

In Willis and Grifﬁths (2008) the authors give the following deﬁnition of the Bistatic Radar Cross
Section.

Deﬁnition 2.8.1 The Bistatic Radar Cross Section (BRCS), denoted by σ

, is a measure of the

energy scattered from the target (aircraft) in the direction of the receiver.

The bistatic radar cross section is a function of bistatic angle β and aspect angle δ and therefore is
more complex than the monostatic Radar Cross Section (RCS).

Considering complex targets we can distinguish three regions of interest for BRCS with respect to
bistatic angle. They are:

i. pseudomonostatic, 0

◦

to 5

◦

ii. bistatic, 5

◦

to 180

◦

iii. forward scatter, ∼180

◦

By the term complex targets we understand a target that is built of discrete scattering centers, such
as ﬂat planes, corner reﬂectors, dihedral with corner not equal to 90

◦

, stationary phase regions.

The ﬁrst – pseudomonostatic region for complex targets is greatly reduced with comparison to
smooth, perfectly conducting targets such as sphere, ogive, cone, for which the region can extend to
β = 100

◦

for a sphere. Therefore the BRCS of complex targets can be reduced to monostatic one,

measured on the bisector of the bistatic angle at a frequency reduced by a factor of cos (β/2), if
the bistatic angle does not exceed 5

◦

. This condition is stated as a variation of equivalence theorem

given in Kell (1965).

2. Used techniques

The second bistatic region starts at the bistatic angle for which the equivalence theorem fails to
determine the BRCS. As previously for complex targets under condition that target aspect angle
is ﬁxed with respect to bistatic bisector, there are three sources of divergence, namely: the phase
between discrete echoing centers now changes; changes in level of energy radiated from discrete
echoing centers; and new echoing centers emerge and some centers disappear.

The third region of BRCS is called the forward scatter and relates to a bistatic angle approaching
180

◦

. Forward scatter conﬁguration was introduced in 1.5.2. In Siegel et al. (1955) it is showed that

when bistatic angle β = 180

◦

the target’s BRCS is relatively small with comparison to the target’s

dimensions and equals σ

= 4πA

/λ

, where σ

denotes BRCS for forward scatter case. In this

case the target might be either smooth or a complex structures.

To obtain the BRCS of some object one must use a suitable technique. Several such techniques have
been developed and among them we can distinguish groups and subgroups as follows:

• Theoretical (Prediction techniques)

– Graphical Electromagnetic Computing (GRECO) (high frequency) (hua QIN and fa WANG,

2002),

– Method of Moments (MoM) (numerical method) (Bhattacharyya, 1991),

– Finite Intergration Technique (FIT) (numerical method) (Weiland, 1996),

– Fast Multipole Method (FMM) (numerical method) (Rokhlin, 1990; Coifman et al.,

1993; Song et al., 1997),

– Shooting and Bouncing Rays (SBR) (Baldauf et al., 1991),

– Closed form Physical Optics (CPO) (Jackson et al., 2010),

– Current Marching Technique (CMT) (Li et al., 2009),

– Finite Element Method (FEM) (Alfonzetti and Borzi, 2000)

• Practical (Measurement techniques)

– Pendulum method (Matsuo et al., 1970)

– outdoor measurements: scaled model, full scale object (Schetne and Mount, 1965; Lane

et al., 1999)

– indoor measurements, most often with use of anechoic chamber: scaled model (Gürel

et al., 2003), full scale object

Arriving at explicit formulae of BRCS for a given object is a very challenging task. Just a couple
discrete shapes have been solved for their explicit formulation. Among them there are sphere and
cylinder.

These shapes were used as an example of explicit formulae of BRCS proﬁles. The ﬁrst one consid-
ered is a perfectly conducting sphere where emitter and transmitter are both polarized perpendicular
to the plane containing the direction of incidence and scattering. A perfectly conducting body in a
time-varying magnetic ﬁeld supports surface currents that shield the magnetic ﬁeld from the interior
of the body. The equations of BRCS for electric and magnetic ﬁelds (σ

and σ

for E-plane and

H-plane, respectively) for the sphere is derived in the following fashion.

2.8 Bistatic Radar Cross Section

(β) =

∞

n=1

(−1)

2n + 1

n (n + 1)

∂P

(cos β)

∂β

+ a

(cos β)

sin β

(2.26)

(β) =

∞

n=1

(−1)

2n + 1

n (n + 1)

(cos β)

sin β

+ a

∂P

(cos β)

∂β

(2.27)

where

))

(2)

)

n−1

) − nJ

)

(2)

n−1

) − nH

(2)

)

(2)

) = J

) − jY

)

(2.28)

, Y

and H

are nth-order Bessel functions (Carrier et al., 2005) of the ﬁrst, second and third

kind, respectively, and P

denotes the associated Legendre function (Filloux et al., 1987) of order

n and degree 1, r

denotes the radius of the sphere and k

, k

2π

is the wave number.

The results of applying (2.26) and (2.27) for different values of the sphere’s circumference k

presented in Fig. 2.9.

In the case of the cylinder the formulations of BRCS for vertical (electric ﬁeld, E-plane) and hori-
zontal (magnetic ﬁeld, H-plane) polarizations are given in (2.29) and (2.30), respectively.

(β) =

2λ

∞

n=0

Ω

(−1)

)

(2)

)

cos nβ

(2.29)

(β) =

2λ

∞

n=0

Ω

(−1)

)

(2)

)

cos nβ

(2.30)

where

Ω

1, n = 0

(2.31a)

2, n = 0,

(2.31b)

the primes denote derivatives with respect to the argument and H

is the Bessel function of the

third kind (Hankel function), r

denotes the radius. To derive the BRCS of a real cylinder the

obtained scattering width needs to be multiplied by the cylinder’s length. It is worth noting that the
BRCS solution is reliable if the wavelength is relatively small with comparison to dimensions of the
cylinder (it also applies to the sphere).

The resulting BRCS of the cylinder for various k

values are illustrated in Fig. 2.10.

In the case of the sphere and the cylinder, inﬁnite sums were approximated by limiting these series
to 2k

and 2k

ﬁrst elements for the sphere and the cylinder, respectively.

2. Used techniques

180

240

300

360

−30

−25

−20

◦

(β

)

= 3

H-plane

E-plane

180

240

300

360

−40

−30

−20

−10

◦

= 6

H-plane

E-plane

180

240

300

360

−40

−30

−20

−10

◦

(β

)

= 12

H-plane

E-plane

180

240

300

360

−40

−20

◦

= 21

H-plane

E-plane

Figure 2.9: Standarized BRCS as a function of bistatic angle of perfectly conducting sphere for
various values of k

2.8 Bistatic Radar Cross Section

180

240

300

360

−30

−20

−10

◦

(β

)

= 19

H-plane

E-plane

180

240

300

360

−40

−20

◦

= 38

H-plane

E-plane

180

240

300

360

−20

◦

(β

)

= 75

H-plane

E-plane

180

240

300

360

−20

◦

= 113

H-plane

E-plane

Figure 2.10: Standarized bistatic RCS as a function of bistatic angle of the cylinder for various
values of k

HAPTER

III

Experiments

This chapter introduces conﬁguration of the Bistatic Radar (BR) system used for acquiring data for
further analyzes.

3.1

Passive bistatic radar scenario

The Passive Bistatic Radar (PBR) scenario used in this work was exploited with one transmitting
site T and one or two receiving sites R. The transmitter used as an illuminator of opportunity, not
synchronized with the receivers, is located in a suburb of Saint Petersburg, Russian Federation. The
receivers were located in the neighborhood of Joensuu, North Karelia, Finland.

Locations of both receiving parties R that were recording Radio Signal Data (RSD), J and M –
ﬁrst letters of the names of the receivers’ operators, and the transmitting station T are depicted in
Fig. 3.2. One of the receiving sites was recording with two antennas, therefore the notation was
given J

and J

, for more details see 3.3. Further we will refer to J as a location rather than a

particular antenna, whereas M will refer to an antenna conﬁguration and location.

Notable distances between receivers J, M – d

, between the transmitter T and receivers J, M –

, d

, are respectively d

J M

= 42.2 km, d

T J

= 301.8 km, d

T M

= 288.4 km.

3.2

Transmitter

The transmitter used for commencing the experiments was the 310 m TV tower built in 1962 with
mounted TV-broadcastning/Frequency Modulation (FM) transmission. The TV transmitter emits on
frequency f

, f

= 49.75 MHz. The Effective Radiated Power (ERP) of St. Peterburg’s TV station

used in this experiment was 149 kW. Until 1985 the tower stood at 316 meters. At that time the
antenna atop the tower was replaced with a more modern version that was six meters shorter.

The transmitting antenna is horizontally polarized and omnidirectional. The pattern of the antenna
is a doughnut so it radiates negligible amount of power in the air space directly above the tower
– the power is uniformly distributed in every direction on the horizontal plane. Because of the
omnidirectional pattern of the transmitter the directional density of radiated power is about 50 %
less than of the same ERP for a dipole (Szóstka, 2006). The omnidirectional antenna is used as a

3. Experiments

dipole is presented in Fig. 3.3. Half dipole length is usually created to tune in to a certain frequency
in which case the length of the dipole is critical. The formulae which describes the relation between
the tuneable frequency, which is referred to as the transmitting frequency f

, and the length of the

dipole L

is presented in (3.1).

150

(3.1)

However the equation for the length of the conductor is only an approximation. Other parameters
that can affect the reception of transmitting frequency are antenna’s height above the ground, nearby
objects and length/diameter ratio of the conductor. To tune in to a transmitting frequency (to achieve
exact resonance), trimming or extending the conductor might be necessary.

75 Ω coaxial

Figure 3.3: Half dipole antenna. I - insulator, R - rope for stability purposes, c

- coaxial splitting

The construction of a 4-element half-wave dipole array is presented in Fig. 3.4. The Yagi used in this
study was located at about 14 m above ground and a gain of 4 to 5 dBd. It is almost omnidirectional
except for the dipole ends which are intentionally directed towards the chosen TV transmitter to
attenuate about 20 dBd of both the TV carrier and the high peak level of signal scattered from
aircraft during the moment of its carrier crossing.

From this point on J

will refer to a receiver with a dipole as the receiving array.

3.3.2

Yagi-Uda

As presented in Straw et al. (2007) the invention of the Yagi-Uda array is attributed to two Japanese
university professors Hidetsugu Yagi and Shintaro Uda from Tohoku University. The Yagi-Uda
array is commonly refered as a Yagi array.

The simplest conﬁguration of a Yagi array consists of two parallel elements: the driven element and
a single parasitic element (Director or Reﬂector). This arrangement is called a 2-element Yagi. If the

3. Experiments

parasitic element is placed opposite to the maximum radiation direction, behind the driven element,
then it is called a reﬂector. Placing the parasitic element in front of the driven element gives it a
director

notation. A Yagi is engineered for Very High Frequency (VHF) and Ultra High Frequency

(UHF) frequencies and can include 30 or more director elements and a single driven element.

The gain and directional patterns of the array are a result of amplitudes and phases of currents
induced in every parasitic element, which on the other hand is given based on the spatial distance
between the driven element and parasitic elements and obviously the length and the diameter of the
elements. This implies that the operational idea of the Yagi array relies on mutual coupling of the
elements.

The following four parameters are relevant factors when it comes to characterizing a Yagi array:

i. Forward gain. Yagi free space gain (over an isotropic radiator in free space) ranges from

5 dBi for a basic 2-element conﬁguration up to 20 dBi for a 31-element UHF design. The
main parameter that determines the gain is speciﬁed in terms of the length of the boom.

ii. Pattern. Gain is strongly determined by the antenna’s directivity pattern, in a particular by

the amount of energy accumulated in particular direction at the expense of energy radiated
in unwanted directions. A 3-element Yagi has a directional gain (main lobe of about 66

◦

for

3 dB loss compared to the peak gain at 0

◦

) of about 7.3 dBi over an isotropic antenna and

5.1 dBd over the half-dipole antenna. Separation of the frontal and rear originated signals is
yet another parameter determined by the directivity pattern. The antenna’s front-to-back ratio
measures the abilities of a Yagi to depress the unwanted rear (90

◦

to 270

◦

) interfering signal.

iii. Drive impedance/Standing Wave Ratio (SWR). The impedance of the driven element is af-

fected by the tuning of the driven element and the spacing and tuning of the parasitic elements.
In some cases inappropriate design of the Yagi can lead to a lack of stability between leading
performance parameters such as SWR, worst case Front-to-Rear ratio (F/R) and gain. In most
cases antennas purposed for high gain usually exhibit large changes of impedance levels with
relatively small changes in frequency. Therefore the SWR changes are wide which may lead
to unwanted changes in the feed cable.

iv. Mechanical strength.

All these parameters must be considered simultaneously to guarantee high receiving performance
and survivability of the array.

The Yagi array presented in Fig. 3.5, was design for reception of the VHF band, particularly for
signal reception of 49.75 MHz. The Yagi antenna will be referred further on as the J

3.4

Radio signal data

Radio Signal Data (RSD) was acquired using antennas described in 3.3.1 and 3.3.2. The receivers
J

, J

and M were tuned to record a spectrum of frequencies that includes the transmitter frequency

, (f

= 49.75 MHz) and accompanying Doppler curves. To obtain the spectrogram form of RSD,

Bracewell (2000), the Short Time Fourier Transform (STFT), presented in 2.2 was used with the
width of a symmetrically positioned Hann window L set to L =∼ 1 s, (8192 samples) and calcu-
lation time step G to G = 0.5 s. With these settings the spectrogram is categorized as overlapped

HAPTER

Long-distance passive multi-static aircraft tracking

In this chapter the multi-static Doppler radar concept based on long-distance Very High Frequency
(VHF) frequency recordings is presented and tested in a two receivers, one transmitter conﬁgu-
ration. The core of the method is a mathematical model of Doppler shift generation in a passive
conﬁguration of many simultaneous transmitting radio stations and receivers that share the Doppler
shift information they have recorded. The tracking capabilities of such a multi-static conﬁguration
are tested in an off-line experiment over a distance of 400 km. Aircraft position was recovered to an
apparent precision of 1500 m. The method also corrects for errors in time synchronization between
receivers.

4.1

Introduction

Air trafﬁc safety is a global concern. Yet air safety radar systems are run to a large extent along
national boundaries. This has caused tragedies at times, like the disappearance of the Air France
ﬂight 447 from Rio de Janeiro, Brazil, to Paris, France, over the Atlantic in 2009, where the avail-
ability of aircraft position over the Atlantic might have resulted in much faster detection of the
debris from the accident. There are examples where long-range aircraft monitoring might even have
prevented accidents, such as the disastrous mid air collisions in Ceritos (CA) 1986, Charkhi Dadri,
India 1996, Überlingen, Germany 2002 and Mato Grosso, Brasil 2006 might have been avoided, if
ground control were provided with correct and exact altitudes and positions of planes.

Part of this safety problem is caused by the inability of monostatic air surveillance radars to see be-
yond the range of their active radar beams. A multi-static radar solution using lower frequencies has
quite some potential in extending the view of air trafﬁc radars to a more global scale. For instance,
Grifﬁths (2010) compares different radar systems and emphasizes in particular the advantages of
Multiple-Input and Multiple-Output (MIMO) radar systems over traditional monostatic ones. At
the moment, the development of transnational multi-static radar systems has been mostly left to
experimental setups. Some recent research devoted to the analysis and testing of multi-static sys-
tems can be found, for instance in Kuschel et al. (2008); Georgescu and Willett (2012) and Doughty
(2008); the issue of the target recognition in the micro–Doppler (µD) radar setup is considered, for
instance in Doughty (2008); Fogle and Rigling (2012).

In this chapter, a novel multi-static air surveillance approach is tested in a multi-static conﬁguration.
It is a passive radar system based on commercial Very High Frequency (VHF) transmitters and

4. Long-distance passive multi-static aircraft tracking

individual radio amateurs. In particular, the DX listener community serves as the receivers. The
DX abbreviation comes from Distance Unknown (DX), and DX’ing is related to activity which is
listening to far away radio stations. The cross section of many pairs of radio stations and many
listeners can create a good global coverage for air trafﬁc monitoring.

We test the viability of such a scenario in an off-line experiment where ﬂight trafﬁc is monitored
over a distance of some 400 km. Two radio amateurs are listening to a remote radio station and
recording the Doppler effect of ﬂights crossing the line of transmission. From a combination of
their shifted Doppler measurements, the position of the aircraft is reconstructed off-line by means
of a mathematical model of the Doppler shift caused by an aircraft that ﬁrst approaches and then
recedes from the line of transmission. The reconstruction of the ﬂight path is compared to the
recordings obtained from the Automatic Dependent Surveillance - Broadcast (ADS-B) transponders
of the aircraft, as reported in the web service http://www.flightradar24.com.

The chapter is organized as follows. The second section explains relation between the paper and pre-
vious work on the subject. The third section presents the mathematical model developed. The fourth
section discusses the capture and processing of the data used in the experiment that is presented in
the ﬁfth section. The chapter concludes with a discussion of the results obtained.

4.2

Previous research on bi and multi-static passive radar

Many authors have contributed methods that are related to the one presented here, but often in a
different context, such as real-time target tracking. In Howland (1999, 1995) the author used a bi-
static radar model with twin Yagi-Uda (3.3.2) antennas together with the Discrete Fourier Transform
(DFT) (2.1) and the difference between resulted phases of received signals in order to ﬁnd estimates
for Doppler and bearings. Target echoes are identiﬁed with use of Cell Averaging – Constant False
Alarm Rate (CA-CFAR) (2.3) and then followed by Kalman Filter (KF) to associate individual
Doppler-Direction of Arrival (DoA) proﬁles. Finally the target’s coordinates and velocity are calcu-
lated with use of Genetic Algorithm (GA), Levenburg–Marquardt algorithm and Extended Kalman
Filter (EKF). The key to approximate target’s track is based on the use of two antennas horizontally
spaced, their phase interferometry is coupled with track estimation process which follows changes
in Doppler shift. Investigation of ground clutter by means of VHF airborne passive Bistatic Radar
(BR) was presented in Brown et al. (2012b). The paper discusses an application of Doppler Beam
Sharpening (DBS) to data collected from airborne receiver which resulted in producing graphical
representation of the area of interest. Successful experiment of detecting air targets by using air-
borne Passive Bistatic Radar (PBR) (1.5.1) with two antennas is presented in Brown et al. (2010).
However in Brown et al. (2012a) the airborne PBR experiment was repeated using a multi-static
radar architecture with two transmitters and compared with ADS-B Mode-S receiver records. The
resulted estimation of air target position was successful too. Ground based long range detection of
up to 700 km of bistatic range with use of Passive Radar Demosntrator (PaRaDe), Frequency Modu-
lation (FM) based PBR, was presented in Malanowski et al. (2012). The work discusses theoretical
limitations of aircraft detection range of PBR taking into account different factors such as Effective
Radiated Power (ERP) of the transmitter, receiver’s dynamic range, size of target, integration time.
Proposition of solving multi-static radar tracking system based solely on Doppler shift information
is demonstrated in Mixon (2012). The work presents strong mathematical basis for solving such a
system but also stresses that precision of measurement of Doppler shifts might prove to be crucial
for performance of presented algorithm.

4.3 Mathematical model

The approach presented here differs from previous study in it that it is genuinely multi-static, with an
arbitrary and possibly variable number of both transmitters - i.e. illuminators of opportunity – and
receivers. The architecture of the system is designed for passive, safety-focused monitoring of air-
trafﬁc by volunteers using peer-to-peer exchange of ﬂight information. Moreover, the detection of
aircraft from VHF Doppler-only signal is conducted with image processing of Doppler images and
mathematical optimization combined. In the current article we describe the algorithm to conduct
such multi-static passive ﬂight trajectory detection and reconstruction and demonstrate its efﬁciency
with a real-life example, albeit an off-line one.

4.3

Mathematical model

This section presents the mathematical model to be applied to the Radio Signal Data (RSD) (3.4) in
order to estimate aircraft trajectories. The model uses spherical coordinates to calculate the distance
between two arbitrary points. All spherical distances are computed along geodesics and derived
from Vincenty’s Inverse Formulae (VIF) (see 2.4) which is accurate enough for the case presented
in this work. However, a technique presented in Karney (2013) can also be used. Hereafter J and M
refer to the receivers and I = {J, M} for notation simpliﬁcation, T stands for the transmitter.

4.3.1

Preprocessing the data

The RSD consists of two time series of signals

(t)

I={J,M}

that have been recorded at two

different locations J and M.

Both sets have the same sample f

rate but a different recording start time t

and different recording

period T

To preprocess the RSD signals the following steps are taken.

First, the signals s

(t) are transformed with the Short Time Fourier Transform (STFT) (Allen and

Rabiner, 1977) into matrix form (spectrogram image) with time-frequency (t, ω) domain expressed
in s and Hz. Namely,

(t, ω) =

N −1

τ =0

[τ ] w

(τ − tG) e

−jωτ

(4.1)

where s

[τ ] is an input signal at time τ , w

(τ ) is a length L window function, S

(t, ω) is DFT of

windowed data centered about time tG and G denotes a hop size (in samples) between successive
DFTs. Hereafter the transformed signals are presented in a form of matrix denoted by S

(t, ω).

After operation presented in 4.1, the frequency domain of matrix S

(t, ω) is shifted as presented in

4.2

ω = ω −

t=1

(t, ω)

= max

t=1

(t, ω)

+ f

(4.2)

4. Long-distance passive multi-static aircraft tracking

where n is a size of matrix S

(t, ω) in time domain, ω denotes new domain of S

(t, ω) and f

the transmitting (carrier) frequency of the transmitted signal. In other words, the column of matrix
S

(t, ω) with the highest average amplitude is subjected to the shift f

. The range of ω remains the

same as the range of ω. This step is important as to ensure that the carrier signal is of the frequency
expressed in f

Next, the matrix S

(t, ω ) is reduced to S

(t, |ω − f

| < f

), where f

denotes the frequency

margin needed to enclose Doppler curves. For the sake of simplicity, we use S

(t, ω

∗

) instead of

(t, |ω − f

| < f

) in the text below.

4.3.2

The model

The theoretical model based mainly on data assimilation comprises the following steps.

• Heuristic approach to spectrograms analyzes based on Doppler history. Both spectrograms are

examined for possible pairs of Doppler curves. Namely, we look for any systematic difference
in time between appearances of Doppler curves on the two spectrograms. The time difference
and distance between receiving stations should coincide with the average cruise speed V

Once the pairs are identiﬁed, we focus on one of them and subject it to a further analysis.

• Time bounds. The region of the spectrograms S

(t, ω

∗

), where the Doppler curve is located,

is cut by time bounds s.t. t ∈

I
l

, t

I
u

. Let n

denote the size of matrix S

(t, ω

∗

) in the

time interval within the bounds t

I
l

, t

I
u

, where t

and t

denote lower and upper time limits

where the Doppler curve was found.

• Checking for swaying of the signal. Since the carrier frequency f

of the signal tends to sway

over time, we need to check its value for the previously deﬁned region. This is done via the
formula in 4.3

∗

t=1

(t, ω

∗

)

= max

∗

t=1

, ω

∗

)

(4.3)

Previously deﬁned value of f

needs to be recalculated due to swaying phenomenon. It is de-

rived by calculating average value of spectrogram S

, ω

∗

) along time domain and searching

for its maximum value. In this way we correct the value of f

. It is important to take this step

into account since further calculations rely on precisely deﬁned value of f

• Edge detection. Next, the matrix S

(t, ω

∗

) is processed with the Canny edge operator (2.5)

with adjusted low and high thresholds. The Canny operator is superior to other edge detecting
methods for this case where the radio signal is fading in and out with time, see Canny (1986).

• Line search within spectrograms. The edge detected image (binary representation of matrix

S) is then used to ﬁnd lines with the Hough Transform (HT). The choice of the HT parameters,
angle θ

and distance ρ

, depends on the shape of the Doppler curve. Finding lines within the

spectrograms is a key element as later they are used as reference lines for synthetic Doppler
ﬁtting. We expect that found segments of lines mostly belong to Doppler curve. In this model,
we assume that Doppler curves emerge from straight-line-trajectories only. This assumption

4.3 Mathematical model

emerges from an idea that system is based on many BR connections thus the network of
bistatic lines is dense enough to ensure good approximation of trajectory even in the case of
aircraft’s change of azimuth. The result of the HT is a family of line segments described in
4.4

t =

sin θ

− ω

∗

cos θ

(4.4)

• The lines detected are used for the ﬁrst estimation of the crossing time. The crossing time

is a time instant at which a family of lines crosses the carrier frequency. The crossing time
is estimated with the secant method and is denoted as t

further. This initial estimation of

the crossing time is not exact because it uses the HT with lines rather than with a model of
the radio signal itself that has a tangent function-like shape. However, it helps to eliminate
unwanted segments that are deﬁned as the ones that satisfy the following condition

s.t.

t < t

I
x

∧ ω

∗

< f

(4.5a)

t > t

I
x

∧ ω

∗

> f

(4.5b)

The segments to be deleted are located on two quarters, ﬁrst and third, of spectrograms with
respect to the origin located at (t

I
x

, f

). This step is responsible for ﬁltering unwanted seg-

ments of lines which presence could have negative effect on a second estimation of crossing
time.

• Extraction of Doppler curves. This step uses results of the HT to extract a Doppler curve from

within the radio signal image S

(t, ω

∗

). First, the family of line segments is smoothed using

weighted linear least squares and a 2nd degree polynomial model that assigns a low weight to
outliers in the regression. Let us denote the smoothed curve as t = f

(ω

∗

) phantomf

. Then,

we look for a set of points from matrix S

(t, ω

∗

) with the maximum value among the points

located close to the smoothed line

∗∗

= ω

∗

max

t, |ω

∗

− f

−1

(t) | < σ

= S

t, |ω

∗

− f

−1

(t) | < σ

(4.6)

where σ

is the standard deviation of differences between the smoothed line f

and the HT

family of lines. Let t = f

(ω

∗∗

)

stand for the newly estimated curve out of this set. The

curve has been found to be discontinuous in most cases of different S

(t, ω

∗

• Removal of outliers. Here, outliers of the curve f

are removed and replaced with interpolated

values. Let us denote the resulting curve as t = f

(ω

∗∗

)

. Then, the second estimation of

the crossing time t

takes place. As before, the secant method is used, but this time f

is the

reference function and t

denotes the new crossing time. At this point, the extraction of the

Doppler curve from matrix S is ﬁnished and the process of ﬁtting starts.

• Search for anchor points. Let us introduce two sets of points, with each set located on a

transmitter T – receivers J, M baselines. The initial location of the sets is dictated by the
minimum and maximum distances that aircraft can traverse during the time span deﬁned by

4. Long-distance passive multi-static aircraft tracking

Establish two sets of points U b

∗

, n

∗

= 1

on transmitter - receivers baselines b

and b

Create trajectories based on geodesics

traversing through pair of points, each from

U b

∗

and U b

∗

. Repeat this step for

every possible pair.

For each trajectory calculate distances d

(t)

and d

(t) deﬁned in 4.8

Based on found distances calculate Doppler

shift f

, f

(t) deﬁned in 4.7

Compare synthetic Doppler curves against the

curve f

Find the pair of anchor points

∗

, ub

∗

that produces the best ﬁt

Is the

accuracy of

ﬁt high

enough?

stop

Produce new set of anchor points based on 4.9

yes

Figure 4.1: Organigram of procedure presented in ’search for anchor points’ step.

J
x

− t

M
x

. The minimum and maximum distance values correspond to the lower and upper

limits that we set on the cruise speed of aircraft V

. As a result, we have two or four solutions,

depending on the time span and distances between the transmitter and receivers. The solution
is then reduced to the pair with the shortest spatial distance to the receivers.

Let us denote these two sets of points as U b

I
a

∗

, n

∗

= 1

and refer to them as anchor

points

. The anchor points are associated with the crossing time so that at the moment of the

crossing time, t

, the aircraft is intersecting the anchor points (one from each set). This

assumption emerges directly from the following equation for f

= 0

(t) =

d (d

(t) + d

(t))

(4.7)

In (4.7), f

(t) denotes Doppler frequency as the function of time, c is the velocity of prop-

agation of electromagnetic waves (light), d

(t) and d

(t), I = {J, M} are distances from

the transmitter to the aircraft and from the aircraft to the receiver, respectively. The distances

4.3 Mathematical model

are derived from the following formulae

TA(AI)

(t) = 2R

+ alt (t)) 1 − cos ξ

TA(AI)

+ alt (t)

(4.8)

alt

Figure 4.2: Diagram that explains geometry used in 4.8

where R

= 6371 km

is the mean radius of the Earth, ξ

and ξ

correspond to great

circle arcs, measured in degrees and connecting the transmitter with the aircraft and the air-
craft with the receiver, respectively. alt denotes the altitude of the aircraft. Here, we assume
that both the transmitter T and the receiver I are within the visible horizon with respect to the
aircraft.

For each point from set U b

J
a

∗

we take every point from set U b

M
a

∗

and create trajectories

based on geodesics traversing through these points. These trajectories are then used to esti-
mate a family of synthetic Doppler curves by using (4.7) and compared against the curve f

The pair of two anchor points

j,1
a,1

, ub

k,2
a,1

that produce the best ﬁt is then used as reference

points for creating new sets of anchor points, as presented in 4.9

a,n











1,1
a

∗

= ub

max(j

−1,1),1

∗

−1

(4.9a)

g,1
a

∗

= ub

min(j

+1,g),1

∗

−1

(4.9b)

1,2
a

∗

= ub

max(j

−1,1),2

∗

−1

(4.9c)

g,2
a

∗

= ub

min(j

+1,g),2

∗

−1

(4.9d)

To clarify, in 4.9 we establish a new sets of points with location between two closest neigh-
bours of previously found points ub

∗

, ub

∗

. This step is repeated recursively, until the

4. Long-distance passive multi-static aircraft tracking

desired degree of precision is achieved. Let us denote the Doppler curve that corresponds
to the best ﬁt after n

∗

= N

∗

iterations by f

and the estimated two sets of anchor points by

U b

I
a

∗

, n

∗

= N

∗

• Extrapolation. In order to ﬁnd a trajectory for the region beyond the anchor points, we extrap-

olate the estimated trajectory using its characteristics, such as velocity and azimuth. Limits
for extrapolation are dictated by the previously deﬁned time limits (t

I
l

, t

I
u

) and the velocity of

the aircraft.

4.4

Processing of an Experimental Data Set

In this section, we describe the data used to demonstrate the application of the mathematical model
presented above.

RSD was recorded by two receivers referred to by J

and M. In addition to RSD, we also retrieved

data from the web service http://www.flightradar24.com Flight Radar 24 (FR24) for
validation of model results.

The recording session for both sets of data (RSD and FR24) took place on July 12, 2012. The
spatial region that was taken to be the scope during FR24 recording encircles the southeastern
corner of Finland and the adjoining part of the Russian Federation. Data was retrieved from the
portal with a sampling rate of 10 s. However, the response from the portal was not immediate and
the sampling time thus varies. RSD was acquired with sampling rate of 2 MS/s. The receivers were
tuned to record spectrum of frequencies that includes transmitter frequency

= 49.75 MHz)

and accompanying Doppler curves. Synchronisation with TV station was not required.

Locations of both receiving parties that were recording RSD J and M and the transmitting station T
are depicted in Fig. 4.5. The ERP of St. Peterburg’s TV station used in this experiment was 149 kW.
The antenna patterns that were used during the recording session include:

• for receiver J rotatable horizontal 4-element 50 MHz Yagi, gain about 5 dBd, height about 7 m

above ground, directive pattern typical to 5-element Yagis, and

• long wire of 250 m, which is slightly directional to northeast and has almost similar back lobe

pattern to south west, for receiver M.

The two baselines that join the positions of the receivers J and M with the transmitting station T are
denoted by b

and b

, respectively, and plotted as dashed lines. The map in Fig. 4.1 also displays

ﬂight trajectories reconstructed from FR24 data. The ﬂights are numbered with respect to the time
instances when they were discovered by the system within the deﬁned region indicated in Fig. 4.1
by a dashed circle. Notable distances between receivers J , M , between receivers and the transmitter
J , S and M , S are respectively d

J,M

= 42.2 km, d

J,S

= 301.8 km, d

M,S

= 288.4 km.

A total of twenty four FR24 ﬂights were recorded. From this set, three ﬂights were classiﬁed to
group G1, namely, ﬂights 3,10 and 11. Group G1 gathers ﬂights which might have left Doppler
traces within RSD signals. We can notice that all the three trajectories end prematurely with respect
to the observed region. This is due to temporary lack of information on aircraft position from ADS-
B located in Joensuu, the eastern-most observation point on the ﬂight corridors used by the ﬂights.

4.4 Processing of an Experimental Data Set

Figure 4.3: Map with recorded FR24 ﬂight trajectories (solid black), locations of listeners J and M
and transmitter S together with two baselines b

and b

between the listeners and the transmitter.

RSD Doppler curve indicated with 1 in Fig. 4.4 that corresponds to FR24 trajectory 10 is taken as
an example of application of the system.

In this work we do not use any prior information on aircraft’s azimuth but we assume that an altitude
is constant for aircraft moving with cruising velocity and is equal to 11 000 m.

Fig. 4.4 presents spectrograms in the form of S

(t, ω

∗

) with frequency margin f

, f

= 199.22 Hz.

Here, on both RSD spectrograms, J and M, we can notice Doppler traces. The Doppler traces are
classiﬁed to be caused by the movement of three different aircraft. The approximate differences in
time between each pair of Doppler curves crossing the carrier frequency f

is equal to 2 min 14 s,

2 min 26 s and 2 min 27 s for Doppler curves 1 to 3 respectively. The Doppler curve crossing the
carrier frequency corresponds to the aircraft crossing baselines between the transmitter and the
receiver. Therefore, the deviation in time differences might come from different crossing locations
or different aircraft cruise speeds.

4. Long-distance passive multi-static aircraft tracking

-199.22

49750000

+199.22

19:52:05

20:27:20

20:47:51

21:27:03

21:38:08

-50

-39.9

-29.7

-19.6

-9.5

-199.22

49750000

+199.22

19:52:05

20:00:04

20:29:34

20:50:17

21:29:30

21:38:08

-70

-58.1

-46.3

-34.4

-22.6

Figure 4.4: Spectrograms for signals s

(t)

I={J,M }

together with marked three pairs of classiﬁed

Doppler curves. Horizontal and vertical axes are in Hz and wall clock units, respectively.

4.5

Case study

Among the pairs of Doppler curves presented in Fig. 4.4 we choose pair number 1 and consider it
as an example for the application of the mathematical model presented in the earlier section. The
next step involves setting lower and upper time limits t

, t

for the HT, as depicted in Fig. 4.5. The

lower and upper limits are set so that the visible RSD Doppler curve lies between them. Then the
Canny edge operator is applied to the limited part of the spectrograms and as a result we have a
binary representation (black and white) of the matrices S

(t, ω

∗

Next, the HT is used to extract segments of lines. Considering Doppler curves related to straight–
line–trajectories and the dimension of a pixel (cell of matrix S

(t, ω

∗

)) which is 0.5 s in the time

domain, or

1 Hz in the frequency domain, the angle parameter θ

is constrained to the interval

(−45

◦

, −1

◦

). The angle parameter should be smaller than −1

◦

to exclude the carrier signal from the

scope of the HT. The result of the HT is presented in Fig. 4.5 as segments of lines, together with the
ﬁrst approximate crossing time t

Next, the segments that fulﬁl the conditions set by 4.5 are removed. After this, HT results are applied
for extracting Doppler curves from the spectrograms; this step is explained in detail in the descrip-
tion of the mathematical model. Extracted Doppler curves are depicted in Fig. 4.6, together with

4.5 Case study

-199.22

49750000

+199.22

20:16:53

20:18:16

20:26:56

20:32:53

20:34:00

-50

-39.9

-29.7

-19.6

-9.5

-199.22

49750000

+199.22

20:16:53

20:26:13

20:29:18

20:32:30

20:34:00

-70

-58.1

-46.3

-34.4

-22.6

Figure 4.5: Results of the HT. The segments are depicted with red colour, time limits t

, t

– with

white dashed lines, crossing time t

– with white thick dashed lines.

the reﬁned estimates of crossing time t

. The differences between the ﬁrst and second estimates

of crossing times are 6 s and 2 s for J and M, respectively. In that regard, the second approximation
might prove crucial for accurate aircraft trajectory estimation.

Fig. 4.7 presents results of the preliminary ﬁtting. The trajectory that corresponds to the Doppler
curve for which the cost function f

− f

has the lowest value is presented on the map. The map also

includes marks for the anchor points (10 points). The points are located so that the full length of the
baselines and a distance of approximately 20 km beyond that is taken into account during the ﬁtting
process. The spectrograms present the ﬁt with cost function values equal to 0.314 Hz and 0.28 Hz
for J and M, respectively. As stated before in the mathematical model, the preliminary ﬁtting is
based on ﬁnding the best anchor points. These anchor points correspond to crossing times in the
time domain. Therefore, the curves visible on spectrograms in Fig. 4.7 are limited by crossing times
t

The preliminary ﬁtting deﬁnes ﬂight parameters, such as velocity, azimuth and location, with respect
to time. With this knowledge, we can extrapolate the trajectory beyond the baselines to meet the
condition t ∈ min t

I
l

, max t

I
u

. After this operation we once again calculate Doppler curves and

ﬁnd the cost function which in this case equals 0.561 Hz and 0.635 Hz for J and M, respectively.

In Fig. 4.8, the extrapolated trajectory and corresponding Doppler curves are presented. The map

4. Long-distance passive multi-static aircraft tracking

-199.22

49750000

+199.22

20:16:53

20:18:31

20:27:02

20:32:58

20:34:00

-50

-39.9

-29.7

-19.6

-9.5

-199.22

49750000

+199.22

20:16:53

20:26:37

20:29:17

20:32:35

20:34:00

-70

-58.1

-46.3

-34.4

-22.6

Figure 4.6: Spectrograms with extracted RSD Doppler curves f

depicted with red solid lines.

Reﬁne crossing times marked with white thick dashed lines.

also includes aircraft trajectories from FR24 for comparison. This trajectory is selected based on a
ﬁne alignment between a synthetic Doppler curve and the estimated Doppler curve f

. Moreover, it

is shown how the Doppler curve f

prolongs the FR24 Doppler curve. However, we can notice that

the latter curve is quite noisy, when compared to the ﬁrst one. Extrapolation of the FR24 trajectory
and its intersection with the two baselines b

results in deriving the location of two points l

. These

points are referred as to anticipated anchor points later in this section.

The fact that two trajectories, namely FR24 and estimated one, overlap over some portion of time is
used to show the spatial proximity between them. Fig. 4.9 presents results of comparing positions
of the two trajectories with respect to time. The location of the aircraft dictated by the estimated
trajectory slightly differs from the FR24. Moreover, the FR24 trajectory seems to wander around
the estimated trajectory. The ﬂight speed is constant in the case of the estimated ﬂight track which
is not the case with the FR24 record, probably because of lack of synchronizity between the various
clocks involved. The difference between their azimuths is 0.22

◦

The results demonstrate that the difference between positions of aircraft reconstructed with both
data sets vary between 200 m and 1600 m along the common ﬂight trajectory. The difference across
the trajectory is less than 100 m.

Further analysis is conducted to ensure that the two sets of RSD data are properly synchronized

4.5 Case study

-199.22

49750000

+199.22

20:16:53

20:18:31

20:27:02

20:32:58

20:34:00

0.31401

-50

-39.9

-29.7

-19.6

-9.5

-199.22

49750000

+199.22

20:16:53

20:26:37

20:29:17

20:32:35

20:34:00

0.28138

-70

-58.1

-46.3

-34.4

-22.6

Figure 4.7: Results of preliminary ﬁtting. Map: found trajectory with red solid line, anchor points
with small black marks. Spectrograms depict Doppler curves corresponding to the trajectory with
red colour. Yellow dashed lines depict previously extracted RSD Doppler curves f

with each other. This analysis shows that in the case of un-synchronized PC clocks between the two
listeners, which results in faulty relative start times t

, we are still able to determine the time shift

needed to synchronize the data. To detect the right time shift we perform a number of trajectory

4. Long-distance passive multi-static aircraft tracking

-199.22

49750000

+199.22

20:10:53

20:18:31

20:27:02

20:32:58

20:34:00

0.561

-50

-39.9

-29.7

-19.6

-9.5

-199.22

49750000

+199.22

20:10:53

20:26:37

20:29:17

20:32:35

20:34:00

0.635

-70

-58.1

-46.3

-34.4

-22.6

Figure 4.8: Result of extrapolation: trajectory and corresponding Doppler curves – red colour. Blue
noisy curves on bottom pictures are Doppler curves resulted from recorded FR24 location/time data.

estimation runs with varying start time t

for listener J. During these runs we check the cost

function f

− f

and the distance from the anticipated anchor points l

to RSD-based anchor points

I
a

∗

, n

∗

= N

∗

. Fig. 4.10 indicates that no time shift is needed for the data considered in this case

study.

4.6 Discussion

20:18:43

20:19:03