A Simulation Model for Detection and Tracking Bio Aerosol Clouds using Elastic Lidar

(1)

Institutionen för systemteknik

Department of Electrical Engineering

Examensarbete

A Simulation Model for Detection and Tracking Bio

Aerosol Clouds using Elastic Lidar

Examensarbete utfört i Reglerteknik vid Tekniska högskolan i Linköping

av

Erika Jönsson

LiTH-ISY-EX--08/4207--SE

Linköping 2008

Department of Electrical Engineering Linköpings tekniska högskola Linköpings universitet Linköpings universitet SE-581 83 Linköping, Sweden 581 83 Linköping

(2)

(3)

A Simulation Model for Detection and Tracking Bio

Aerosol Clouds using Elastic Lidar

Examensarbete utfört i Reglerteknik

vid Tekniska högskolan i Linköping

av

Erika Jönsson

LiTH-ISY-EX--08/4207--SE

Handledare: Ove Steinvall

foi, Linköping

Christian Lundquist

isy, Linköpings universitet

Examinator: Thomas Schön

isy, Linköpings universitet

(4)

(5)

Avdelning, Institution

Division, Department

Division of Automatic Control Department of Electrical Engineering Linköpings universitet

SE-581 83 Linköping, Sweden

Datum Date 2008-12-12 Språk Language 2 Svenska/Swedish 2 Engelska/English 2 Rapporttyp Report category 2 Licentiatavhandling 2 Examensarbete 2 C-uppsats 2 D-uppsats 2 Övrig rapport 2

URL för elektronisk version

http://www.control.isy.liu.se http://urn.kb.se/resolve?urn=urn:nbn:se:liu:diva-ZZZZ ISBN — ISRN LiTH-ISY-EX--08/4207--SE

Serietitel och serienummer

Title of series, numbering

ISSN

—

Titel

Title

En Simuleringsmodell för Detektering och Följning av Biologiska Ämnen med Hjälp av Elastisk Lidar

A Simulation Model for Detection and Tracking Bio Aerosol Clouds using Elastic Lidar Författare Author Erika Jönsson Sammanfattning Abstract

Discharges of warfare bio aerosol clouds are powerful weapons in war and terror situations. A discharge of a small amount of a contagious substance can obliterate large areas. The discharges can usually not be seen with bare eyes, hence some tool needs to be used to find bio aerosol cloud discharges. One way is to use a lidar for the detection of clouds. By sending out a laser pulse into the atmosphere some of the light is scattered back. By measuring the backscattered light, the aerosol structure of the atmosphere can be obtained. If a cloud is hit by the laser beam, an increase of light is observed and the cloud can be detected and tracked. In this thesis a tool for simulating the elastic backscattered light has been devel-oped. A graphical user interface for easier handling has also been develdevel-oped. For automatic detection of clouds some detection algorithms have been tested. An-other graphical user interface for presentation of the simulated lidar signals has also been constructed.

The simulator is working well. A lot of different parameters can be changed for the lidar system, the atmosphere and the cloud type. The model works as a helpful tool for specifications of an elastic lidar system to be developed and also as a guide for expected system performance and ways to present the results.

Nyckelord

(6)

(7)

Abstract

Discharges of warfare bio aerosol clouds are powerful weapons in war and terror situations. A discharge of a small amount of a contagious substance can obliterate large areas. The discharges can usually not be seen with bare eyes, hence some tool needs to be used to find bio aerosol cloud discharges. One way is to use a lidar for the detection of clouds. By sending out a laser pulse into the atmosphere some of the light is scattered back. By measuring the backscattered light, the aerosol structure of the atmosphere can be obtained. If a cloud is hit by the laser beam, an increase of light is observed and the cloud can be detected and tracked. In this thesis a tool for simulating the elastic backscattered light has been devel-oped. A graphical user interface for easier handling has also been develdevel-oped. For automatic detection of clouds some detection algorithms have been tested. An-other graphical user interface for presentation of the simulated lidar signals has also been constructed.

The simulator is working well. A lot of different parameters can be changed for the lidar system, the atmosphere and the cloud type. The model works as a helpful tool for specifications of an elastic lidar system to be developed and also as a guide for expected system performance and ways to present the results.

(8)

Sammanfattning

Utsläpp av biologiska ämnen är kraftfulla vapen som används i krigs och terrorsit-uationer och endast en liten mängd smitta kan utplåna stora områden. Utsläppen kan vara svåra att se, därför behövs det hjälpmedel för att hitta utsläppen. Ett sätt är att använda ett så kallat lidarsystem för detektering av moln. Då en laser-puls skickas ut i atmosfären kommer en viss del av ljuset att spridas tillbaka till lidarsystemet. Genom att mäta mängden av det bakåtspridda ljuset kan man få en bild av hur strukturen i atmosfären ser ut. Om ett moln träffas av laserstrålen kommer det bakåtspridda ljuset att öka. Molnet kan då detekteras och följas. I denna uppsats utvecklas ett simuleringsverktyg för simulering av bakåtspritt ljus. För att på ett lätt sätt kunna göra simuleringar har ett grafiskt användargränss-nitt utvecklats. Några olika detekteringsalgoritmer för automatisk upptäckt har utvecklats och testats. Ytterligare ett grafiskt användargränssnitt har utvecklats för presentation av de simulerade lidarsignalerna.

Simuleringen fungerar bra. Flera olika parametrar kan ändras för lidarsystemet, atmosfären och molntyper. Modellen fungerar som en hjälp för att specificera hur ett framtida lidarsystem ska konstrueras och hur resultat ska presenteras.

(9)

Acknowledgments

I would like to thank Ove Steinvall, my supervisior at FOI in Linköping, for let-ting me do this thesis. Thanks to Ove Gustafsson for proofreading of the report, for always having time to read and to give me commens on my writing. I would like to thank Rolf Persson for helping me with all administration during my first days at FOI. Thanks to all people at FOI for making my time pleasant.

Thanks to Christian Lundqvist, my supervisior at ISY for helping me with the report and for reading it and giving me comments. I would also like to thank my examiner Thomas Schön at ISY.

Thanks to my family and friends for listening to me when I needed to talk. Last but not least I would like to thank my beloved fiancé Martin for always beeing by my side and for supporting and encouraging my all way through this work. Thank you!

(10)

(11)

Introduction

This work is made on commission from FOI in Linköping, Sweden. This thesis is divided into two parts. In the first part a simulation tool for simulating the detection of lidar signals hitting bio aerosol clouds in the atmosphere is developed. In the second part some detection algorithms are developed and tested on the simulated lidar signals.

1.1 Problem Description

Bio aerosol warfare agents are extremely dangerous to be exposed to. Hence, a warning system that detects the discharges far away is desirable. Since a long time ago FOI has developed and used atmospheric lidars. Lidar stands for LIght Detection And Ranging. At the moment, one of the existing lidar systems is used for investigation of the aerosols within a European network. Aerosols is a suspen-sion of fine solid or liquid droplets in a gas, which can be seen in the atmosphere. FOI has got a commission to develop a lidar system for bio aerosol detection and classification. Because of this, a simulation tool, detection algorithms and presen-tation are needed. In Figure 1.1 the final system is illustrated. Every laser pulse transmitted scatters back an intensity dependence echo which corresponds to the particle density variation in the beam direction. Those echoes are presented both individually and also like radar, in a PPI (Plan Position Indicator)-image.

1.2 Purpose

The purpose with this master thesis is to

• Develop a simulation model for lidar echoes. The development of the simula-tion model can be based on an existing model which needs improvement with respect to the graphical user interface (GUI) and the possibility to simulate more complicated realeases.

(14)

2 Introduction

• Develop a presentation tool to present the simulated lidar signals in a PPI-image.

• Test and develop a number of detection algorithms for detection of clouds in both single shots and in the whole PPI-image.

• Track a detected cloud using a more narrow angular scanning.

Figure 1.1: An overview of the final system. The laser will sweep back and forth to look for discharges [7].

1.3 Outline

The structure of the thesis is described in this section.

Chapter 2, Lidar Theory presents the lidar equation and all terms in the

equation are explained.

Chapter 3, Simulation of Bio Aerosol Discharges presents all steps in the

simulation process.

Chapter 4, Detection and Tracking of Simulated Bio Aerosol Clouds

describes the detection algorithms and the tracking algorithm.

Chapter 5, Example of Results gives example of results from some

simula-tions and from the detection algorithms.

Chapter 6, Concluding Remarks gives a summary of the thesis and the

results are discussed. Some future directions are given.

Appendix A, Users Manual: Simulation GUI contains some screen shots

(15)

1.3 Outline 3 Appendix B, Users Manual: Presentation GUI contains some screen shots

(16)

(17)

Chapter 2

Lidar Theory

Lidar is commonly used for aerosol and gas detection in the atmosphere. A li-dar system works in a similar way as a rali-dar. The main difference is that lili-dar emits light pulses and radar uses radio frequencies. When light from a laser beam hits molecules and particles in the atmosphere it will be absorbed and scattered in all directions. An amount of the scattered light will come back towards the lidar system and will be collected by a telescope and focused on a photodetector that measures the amount of light backscattered as a function of the range, see Figure 2.1 [1].

Figure 2.1: When a laser hits a cloud the light will be scattered in all directions. Only the light scattered backwards is the light detected by the telescope.

A lidar is able to detect matter that reflects the transmitted electromagnetic energy. In situations when the backscattered light has the same wavelength as

(18)

6 Lidar Theory

the emitted light it is called elastic scatter. This is the scattering used in this thesis [2].

A pulsed laser is used from which the range can be received. After a laser pulse has been sent out, some light will be scattered back into the telescope. The detected signal corresponds to a range profile consisting of the backscatter from the atmosphere. To avoid backscattering with high intensity from short ranges the lidar and the receiver are not coaxially configured. Because of this the receiver will receive maximum signal strength some meters away from the lidar system [2].

2.1 The Lidar Equation

The lidar equation is a fundamental equation in laser remote sensing field to relate the received light power with the transmitted laser power, the light transmission in the atmosphere, the physical interaction between light and objects, the lidarsystem efficiency and geometry, etc. One way to express the lidar equation is

P (λ, R) = PL· ηoptical· c · τL 2 · O(R) R2 · βtot(R, λ) · T 2 tot(R, λ) + n(R), (2.1)

where P (λ, R) is the received optical power, PL is the power of the laser, ηoptical

the optical efficiency, c the speed of light, τL the pulse length of the laser, O(R)

the overlap function, R the range, βtot(R, λ) the backscatter coefficient, Ttot2 (R, λ)

the transmission and n(R) the noise [8]. In the next sections all the terms in the equation will be explained more closely.

2.1.1 Power of the Laser P

L

The peak power of the laser is calculated as

PL=

E τL

, (2.2)

where E is the pulse energy of the laser and τL is the pulse length of the laser

light. The pulse length is defined as the full width at half maximum (FWHM).

2.1.2 Optical Efficiency η

optical

The optical efficiency is chosen by the user. This efficiency incorporate an overall efficiency due to losses by optical components in the transmitter and receiver chain.

2.1.3 Range Resolution

c·τL

2

The term c·τLis calculating the length of the volume illuminated by the laser pulse

at a fixed time. The term 1₂ appears because of the light travelling both forth and back into the atmosphere. When the lidar signal is detected at an instant time

(19)

2.1 The Lidar Equation 7

t after the leading edge of the pulse it comes from the distance R1 = c · t/2.

At the same time, light produced by the trailing edge arrives from the distance

R2 = c · t/2. This leads to ∆R = R1− R2 = c · τL/2 which is the length of the

volume from which backscattered light is received at an instant time and this is called the effictive pulse length [8].

2.1.4 The Overlap Function O(R)

In the lidar system the laser and the telescope are biaxial, which means they are parallel as can be seen in Figure 2.2. Because of this a function called the overlap function needs to be taken into account. The received lidar signal is attenuated at short ranges due to insufficient overlap between the laser lobe and the receiver field of view [1].

Figure 2.2: The laser and the telescope are parallel, which leads to a variation in the overlap between the field of view for the laser and the telescope.

Because of this, the laser cannot be completely imaged onto the photodetector for short ranges since only a part of the actual lidar signal is measured. The overlap varies with the distance and depends on the laser beam diameter, the shape and divergence, the telescopes imaging properties, the receiver field of view, and the location of the emitter and receiver optical axes relative to each other. The overlap function, which can be seen in Figure 2.3, starts with the value zero at the lidar and is thereafter becoming a unity when the laser beam is completely overlapping with the telescope field of view [8].

(20)

8 Lidar Theory

Figure 2.3: The overlap, O(R), between the laser and the telescope, as a function of the range, R [8].

2.1.5 Backscatter Coefficient β

tot

(R, λ)

The backscatter coefficient, βtot, is a measure of the atmospheric backscatter into

the receiver. For the elastic scattering the wavelength of the returning photons are conserved, only the direction of propagation is modified after interaction with aerosols and particles in the atmosphere. When the light hits a particle or an aerosol in the atmosphere scattering occurs at several angles, this scattering de-pends upon the geometric size and shape of the particle, the refractive index of the particle, the wavelength of the incident light, and on the particulate number density [1].

The backscatter coefficient can be expressed as

β(λL, λ, R) = N (R)

dσs_(λ L, λ)

dΩ , (2.3)

where N (R) is representing the concentration of the substance with the differential backscatter cross section dσs_(λ

L, λ)/dΩ. The crossection, σs, represents the cross

sections scattering at the wavelength λ when the incident wavelength is equal to λL

which is equal to the laser wavelength [2]. The amount of light scattered forward is greater than the light scattered backwards, here it is only the light scattered right backwards, at 180◦, that is of interest in lidar measurements, see Figure 2.4 [1].

The backscattering coefficient, βtot, can be divided into two different parts, the

(21)

2.1 The Lidar Equation 9

Figure 2.4: It is only the light scattered at 180◦ that is of interest in lidar mea-surements.

2.1.6 Atmospheric Transmisson T

2

tot

(R, λ)

When the laser beam is going through the atmosphere it is interacting with gases, particles and aerosols in the atmosphere. This interaction is leading to beam de-formations and an increasing divergence of the laser beam. After every interaction a fraction of light is lost due to reflection and absorption and the beam will lose its intensity. The extinction is mainly depending on the concentration of the particles and on the particle size in the atmosphere [2].

The transmission is a function of the range and is calculated for both the transmitted and the scattered light as

T_tot2 (R, λ) = T (λL, R)T (λ, R) = exp[−

R

Z

0

αtot(λL, R) + αtot(λ, R)dR]. (2.4)

The transmission factor of the atmosphere can be described with the extinction coefficients αtot(λL, R) and αtot(λ, R), where λLis the wavelength for the outgoing

light, λ is the wavelength for the incomming light and R is the range. When elastic scattering occurs the wavelength for the outgoing and incoming light are the same, i.e. λL = λ [2]. The extinction coefficient is the sum of all the transmission

losses and can, as the backscatter coefficients, be divided into two different parts,

αtot= αcloud+ αatm [2].

2.1.7 Noise n(R)

The detected lidar signals will always have some amount of disturbing noise. In daytime the background noise is dominated by direct or scattered sunlight, and at night the moon and stars together with artificial light sources will disturb the signal received. Even the detector is a source where undesired noise is produced, in Section 3.1.4 the simulation of the noise is described [8].

(22)

(23)

Chapter 3

Simulation of Bio Aerosol

Discharges

In this chapter the simulation model will be described and presented. For the simulation the lidar equation, discussed in the previous chapter, is used and in this chapter all steps in the simulation process will be explaind. For easier simulations a simulation GUI has been developed, a users manual and screenshots can be found in Appendix A.

3.1 The Simulation Model

When simulating the lidar signal, the lidar equation is used for calculating the received power. From Equation (2.1) the separate terms can be simulated.

3.1.1 The Overlap Function

The overlap function looks like the example given in Figure 2.3. In the simulation, a similar function is constructed. From the start the function is equal to zero. The length of the zero value is determined by the user, a default value for the blocking range is 150 m. Another value for the user to choose is when the overlap function is going to reach 50% overlap, the default value is set to 260 m.

3.1.2 Backscatter

The backscatter coefficient, βtot, consists of backscatter from the atmosphere,

βatm, and from the cloud, βcloud as said before. When a wavelength of 1.5 µm is

used the backscatter from the atmosphere is approximated by

βatm = (23/V ) · 7.2 · 10−7 [msr−1], (3.1)

(24)

12 Simulation of Bio Aerosol Discharges

as

βcloud= βc· Ncloud [msr−1], (3.2)

where Ncloud is the cloud concentration in particles/m3, see Section 3.3, and βc

is a constant depending on the radius and the geometric standard deviations of the particles of the cloud. If the particles have a longer radius and have a high concentration the backscatter coefficient, βc, will be larger compared to the case

with shorter radius and lower concentration. The constant βcis interpolated from

Table C.1 in Appendix C, relevant for the substance BG (Bacillus atrophaeus) [7].

3.1.3 Transmission

The extinction coefficient, αtot, is divided into two different parts, one for the

atmosphere, αatm, and one for the cloud, αcloud. When calculating the extinction

coefficient in the atmosphere, if the laser wavelength of 1.5 µm is used, the equation

αatm = (23/V ) · 0.0382 · 10−3 [m−1] (3.3)

is used, where V is the visibility in kilometres [7]. For the cloud the extinction coefficient, αcloud, is equal to

αcloud= αc· Ncloud [m−1], (3.4)

where Ncloud is the same cloud concentration in particles/m3 as in Section 3.1.2

and αc is interpolated from Table C.2 in Appendix C. The transmission vector is

then simulated in Matlab as

T_tot=exp(-2 * cumsum(alpha_tot)),

where αtot= αcloud+ αatm.

The backscatter, βc, and extinction, αc, coefficients for the cloud can be

deter-mined depending on the radius and the standard deviation of the particles. For calculating the coefficients a table is used where some values are stored depending on radius and standard deviation. By interpolating in the table an appropriate value for the coefficients can be received, see Example 3.1.

(25)

3.1 The Simulation Model 13 Example 3.1

Particle size radius = 0.925 µm

Geometric standard deviaton = 1.3 µm

An approximated value for βc is received by interpolate in the table by using the

Matlab command interp2

beta_c = interp2(X, Y, excel_table, 0.925, 1.3),

where X and Y are a meshgrids and excel_table is the table from Appendix C.

3.1.4 Noise

The noise term is consisting of electronic noise and atmospheric noise.

Electrical Noise

In the detector there is some dark current. This current is flowing in the detector in the absence of any signal or background light. The detector shot noise is generated by random fluctations in the total current. The shot noise current is given by

Ishotnoise=

p

2 · q · IDC· ∆f [A], (3.5)

where IDCis the DC background current due to sun or terrain reflectivity into the

receiver, q is the electric charge of an electron and ∆f is the bandwidth of the pho-todetector [1]. The DC background current is calculated differend depending on if there is sunlight scattering or just scattering from the terrain. The calculations for the IDC due to the sun is calculated as

IDC, sun= U GR ·

π · d2 pix

4 · ∆λ · Edet [A], (3.6) where UGR is the unity gain responsitivity, Edetis the irradiance in the detector

plane, dpix is the diameter of the photodetector and ∆λ is the bandwith of the

optical filter for the photodetector.

For calculation of the IDC due to terrain scattering the following equation is

used IDC, terrain= U GR · π · d2 pix 4 · ∆λ · Esun· Topt· π · D2 4 · ρterrain π [A], (3.7)

where Esun is the sun spectral irradiance, Topt is the optical efficiency, D is the

diameter of the telescope, and ρterrain is the terrain reflectivity.

(26)

14 Simulation of Bio Aerosol Discharges

Idetnoise=

p

2 · q · ID· ∆f [A], (3.8)

where ID= 1 · 10−9 is the unity gain dark current for the detector.

The total noise is then given by

n = p(I

2

shotnoise+ Idetnoise2 ) · M2· Enoise+ IP N2

U GR · M · Topt· ηe

, (3.9)

where M is the detector gain, Enoisethe excess noise and IPNis the preamplitude

noise current and ηe is the efficiency factor.

Atmospheric Noise

By studying a number of real data a way of simulating the natural variations in the atmosphere has been developed. The variations are approximatly Gaussian distributed. To simulate the variations the Matlab command randn is used. A randomly distributed vector with zero mean and a chosen standard deviation is constructed and smoothed with a moving average as

n_atmosphere=smooth(std*randn(size(R)), SPAN, moving),

where the standard deviation and SPAN are chosen by the user. The SPAN determining the size of the groups when the moving average is applied and the standard deviation is given in percent. The backscatter from the atmosphere,

βatm, is then modelled with the atmosphere vector as

βatm = βatm· (1 + n_atmosphere). (3.10)

3.2 Calculation of the Range in Clear Atmosphere

In clear atmosphere the range sensitivity is different depending on the system SNR and on the visibility. By using the lidar equation and solve it for the range an ideal range can be calculated. Cloud is easier to detected within the calculated range then beyond the range. Beyond the calculated range the cloud will be mixed with the atmospheric noise and it is harder to discover. The range is calculated with the Lambert W function as

Rv= W (αatm· √ S) αatm , (3.11)

where Rv is the range, W is the Lambert W function and S is equal to

S =ηoptical· PL· βtot· c · τL· Ar

(27)

3.3 Cloud Concentration 15

where Ar is the receiver area and SNR is the signal to noise ratio for the

atmo-sphere [6].

3.3 Cloud Concentration

When simulating the cloud concentration a model made by Leif Persson at FOI in Umeå has been used [3]. The model simulates different discharges of bio aerosol clouds. In one case there is a continuous plume discharge and in another case there is a puff discharge. In the Umeå model the concentration from the cloud is received as Ncloud(x, y, z) = Q U · (2π)−3/2_σ xσyσz · (exp(−x 2 2σx − y 2 2σy − z 2 1 2σz ) + (exp(− z 2 2 2σz )),

where σxand σy are the standard deviations for the cloud radius in the horizontal

and vertical directions in meter. Q is the number of particles per second and U is the wind velocity in meters per second. The standard deviations are given by the analytical formulas of the form

σ(x) = α(x + x0)

(1 + β(x + x0))γ

, (3.13)

where the constants α, β and γ are given by analytical formulas, depending on Pasquill stability class, averaging time and ground roughness. For more informa-tion about this equainforma-tion, see [3].

(28)

(29)

Chapter 4

Detection and Tracking of

Simulated Bio Aerosol

clouds

In this part algorithms for automatic detection of clouds are presented. The simu-lated data presented in Chapter 3 is used for evaluation of the algorithms. A GUI for presentation has been developed, a users manual with instructions and screen shots can be found in Appendix B.

4.1 Pre Processing

Before detection and tracking algorithms can be applied there are some pre-processing steps to go through.

4.1.1 Range Correction

The lidar equation (2.1) is divided by R2_{, because of this the received signal}

needs to be multiliped with R2 _{to compensate for the range dependency. In the}

measurement a range resolution of 3.75 m is used, which means there is a range of 3.75 m between each sample. Therefore every shot is multiplied with a range vector to compare for the range dependency.

Raw data from one shot can be seen in Figure 4.1. In the figure two local maxima can be recognized, one at 911 m and the other at 1039 m. These peaks are detected clouds from a discharge at a field trial in Umeå 2008-09-24. In the beginning of the signal the typical appearance can be seen and is caused by the overlap and the transmission function. In Figure 4.2 a range corrected signal can be seen.

(30)

18 Detection and Tracking of Simulated Bio Aerosol clouds

Figure 4.1: Raw data from one single shot. Received in Umeå in September 2008.

4.1.2 Moving Average

A simple way to smooth noisy data is to use a method called moving average. In this method a set of data points is analysed by taking the average of groups of the data. This is a commonly used method to smooth out short term fluctuations from noisy data sets. Depending on the size of the average groups more or less of the fluctuations will be removed. The moving average is calculated as

st= 1 k k−1 X n=0 xt−n= xt+ xt−1+ xt−2+ · · · + xt−k+1 k = st−1+ xt− xt−k k ,

where stis the mean of the last k observations. A smoothed signal can be seen in

Figure 4.3.

4.2 Change Detection

By sweeping the laser sideways between the shots, typically at an angular stop of one or a few milliradians, a PPI-image can be constructed. When a number of images from the same scene has been received some image processing techniques can be used for finding clouds. One technique to find changes between two images is called change detection. The algorithm will find anomaly behaviour by comparing two images to each other and if changes occur an alarm will be activated. Changes

(31)

4.2 Change Detection 19

Figure 4.2: A range corrected single shot. Received in Umeå in September 2008.

in images can be caused by appearance and disappearance of objects, movement of objects relative to the background or shape and changes of objects [5].

A simple way to use change detection is to take two images P1 and P2 and

calculate the difference in intensity between the two images. Every pixel/sample in the image represents the intensity at that specific position. The images P1 and

P2 must consist of two sweeps with the laser at the same position and the shots

made in the same angles for both sweeps. The difference image will be calculated as

D(xi,j) = P2(xi,j) − P1(xi,j), (4.1)

where xi,j represents the pixel coordinate, P1(xi,j) and P2(xi,j) represents the

intensity value at location xi,j in the two images P1 and P2 and i = 1, 2, . . . , n

and j = 1, 2, . . . , n represents the index of the pixels. A new image D(x) will be received in which the changes, if any, will be represented. In the differenced image a threshold must be set to determine if any noticeable changes has occured. Due to the natural variations in the atmosphere few or none of the differences will have values equal to zero, hence a threshold is needed. All the values in the difference image are compared to a threshold T1 and values above the threshold are set to

be active, i.e. A(xi,j) = 1 if D(xi,j) > T1 0 if D(xi,j) < T1. (4.2)

(32)

Figure 4.3: Moving average applied on a range corrected signal. The upper image shows the original signal and the lower image shows the averaged signal.

Thresholds are always tricky to determine. In this work the threshold T1 is

de-pendent on the mean of the differences image. The mean is calculated and the threshold is set to be a chosen value times larger than the mean value.

In the next step the number of active pixels will be counted and compared to another threshold T2. If the numbers of active pixels are reaching the threshold,

then the image can be considered to be an active image; otherwise the image is inactive, meaning no noticeable changes have occurred. If the size of the image is

n x m the sum of the active pixels can be calculated as

λ = n X i=1 m X j=1 A(xi,j). (4.3)

To determine if the image is an active image the sum of the pixels λ is compared to the threshold T2 as

Active =

1 if λ ≥ T2

0 if λ < T2.

The threshold T2will be designated as a percentage of changed pixels compared to

the number of pixels in the whole image [9]. This tool will be used when searching for interesting targets to examine closer.

(33)

4.3 Single Waveform Analysis 21

Algorithm 1 Change Detection

1. Calculate the difference between two following sweeps.

2. Compare to the threshold T1 and set values above the threshold to one and

values below the threshold to zero. 3. Calculate the active pixels.

4. If the sum of the active pixels are above the threshold T2 set an alarm flag

to the value one.

4.3 Single Waveform Analysis

In this part single shots will be used for cloud detection. By examining every single shot local maxima can be found and stored in a matrix. Thereafter all shots in the sweep can be taken into account and local maxima close to each other are supposed to belong to the same cloud.

4.3.1 Finding Maxima

Single shots will be examined separately in this phase. A typical range corrected shot will look like Figure 4.2. In this shot two hits are registered, as said before. To find the hits, e.i. the local maxima, in the signal a peak detector will be used.

4.3.2 Simple Peak Detector

To detect peaks in a signal a simple peak detector using a fixed threshold can be used. A problem with this method is how to find the optimal threshold. All samples in the signal will be compared to the threshold and if a sample is exceeding the threshold an additional condition must be fulfilled. The exceeded sample has to be a local maximum to fulfil the quality to be labelled as a peak. This condition is fulfilled if the next sample, x[n + 1], and the previous sample, x[n − 1], have a lower value than the current sample, x[n].

The algorithm for the simple peak detector is illustraded in Figure 4.4. All peaks found in the signal are stored in a matrix where information of sample numbers and values are stored. When a sequence of shots are examined the shot number will be stored in the matrix as well.

In Figure 4.5 the local maxima beyond a threshold are marked with stars and the threshold used is marked with a dotted line. Before the peak detection the signal was range corrected and filtered with a moving average. The maxima detected are located at (911.3m, 0.8213) and (1039m, 2.208).

The threshold can be set as a certain percentage above the standard deviation of the noise. To determine the threshold some empty shots are needed. A mean value of about ten shots is used and the standard deviation and mean value of the

(34)

Figure 4.4: A flowchart of the peak detection process.

noise are calculated. Only samples from a range of about 200 m are used. An example of a threshold calculation can be seen in Figure 4.6.

4.3.3 Removing False Peaks

In the signal there can be interferences, like the laser hitting a bird or some random high peaks which are not due to a cloud hit; all of them are false hits and are not wanted. If those hits reach the threshold for being detected those tops are detected as peaks, which is not desirable. The false peaks cannot be found only by looking at one single shot. To know were the false peaks are situated a number of shots must be compared; after that some of the false peaks can be reduced from the stored peak matrix. If the found maxima are located at about the same range in contiguous shots they are probably belonging to the same cloud.

In Figure 4.7 there are some outliers detected by the peak detector. There are some peaks which do not belong to a cloud. To get rid of the outlier the euclidian distance between all peaks are determined. The Euclidian distance between two points P = (px, py) and Q = (qx, qy) is calculated as

D =k P − Q k=

q

(px− qx)2+ (py− qy)2, (4.4)

it is the distance in samples between the points. If the distance, D, between one peak and all the others are longer than a predetermined length, the peak is taken away from the stored matrix, see Figure 4.8.

4.3.4 Finding the Left and Right Edge of the Cloud

When the false peaks have been removed the left and right edge of the cloud needs to be found. This is important information when tracking a cloud. To determine the left and right edge of the cloud a shortest allowed euclidian distance, see Equation (4.4), between peaks belonging to the same cloud is determined. Therafter the first peak from the left in the image is used, and the euclidian distance to the next peak in the image is calculated. If the distance is shorter than the predetermined distance the peaks belong to the same cloud. If the distance

(35)

4.4 Finding the Centre of Mass for the Cloud 23

Figure 4.5: Local maxima for a single shot. The maxima is marked with stars and the threshold used is marked with a dotted line.

is too long the peaks are beloning to different clouds. By going through all peaks the left and right edges can be found.

Algorithm 2 Finding the Left and Right Edge of the Cloud

1. Find local maxima for every single shot, see Section 4.3.2. 2. Remove outliers, see Section 4.3.3.

3. Calculate the Euclidian distance between all local maxima and determine the right and left edge of the cloud, see Section 4.3.4.

4.4 Finding the Centre of Mass for the Cloud

In this section the whole sweep is used to calculate the centre of mass for the simulated cloud.

4.4.1 Binary Image

A binary image is an image consisting of only two possible values for each pixel. Typically the two different values used are 0 and 1, which correspond to black and

(36)

Figure 4.6: The threshold is determined as a certain percentage above the standard deviation of the noise.

white. The zeroes, the black area in the image, are used to imagine the background and the ones, the white areas, are used to picture the foreground. This method can be used to show important things in an image, to show the pixels with values over a certain threshold.

By determining a threshold as described above the image of one sweep can be transformed into a binary image. After this operation has been run the areas with higher intensity are standing out from the background as white objects.

In the binary image there might be noise together with the important struc-tures, see Figure 4.9. As a human it is easy to separate the large white area from the small ones caused by the noise, but it is harder to implement this function in software. One way to separate the cloud from the noise is to use a morpho-logical labelling. A Matlab function called bwlabel is connecting components in a binary image. Every white area is labelled with a number depending on its neighbours. Now it is possible to access all white areas separate from each other. The biggest area of ones is supposed to be the cloud area.

By filtering the image with the received mask only the desired cloud is left in the scene. Now the centre of mass can be calculated, this is described in the next section.

4.4.2 Centre of Mass of the Cloud

In a rigid body the centre of mass is a fixed point in the object. The point must not lie in the object itself, but it is fixed in the same place relative the body. For

(37)

4.5 Tracking the Cloud 25

Figure 4.7: The peaks detected are marked with stars. The right image is an enlargement of the detected cloud.

a cloud, or other variable structure, the centre of mass is moving when the body is moving and changing its structure.

The definition of the centre of mass is defined as

x = PN i=1mixi PN i=1mi , y = PN i=1miyi PN i=1mi ,

where m is the mass of the object [4].

4.5 Tracking the Cloud

When a cloud is detected the system is supposed to make the sweep narrower. If the sweep is narrower the laser can sweep faster and the precision will be improved. The shots can be closer to each other which will improve the resolution.

(38)

Algorithm 3 Track

1. The function Change Detection, see Algorithm 1, is used to compare two sweeps with each other. If any changes has occured a flag, Target Hit will be set to the a value of one.

2. If a target hit occurred the left and right edge of the cloud will be found by using the function Find Edges, see Algorithm 2. If no changes occurred go back to 1.

3. The sweep will be narrower using the left and right edge of the cloud. 4. Go back to 1.

Figure 4.8: Some of the false peaks have been removed, in compresion with Fig-ure 4.7

(39)

4.5 Tracking the Cloud 27

Figure 4.9: An example of a binary image with a cloud in the field of view. Same sweep as shown in Figure 4.7 and Figure 4.8

(40)

(41)

Chapter 5

Example of Results

In this chapter example results from some simulations will be presented. All simulations are run with the default values, found in Figure A.2, if nothing else is written.

In Figures 5.1, 5.2, 5.3 and 5.4 single shots in an empty atmosphere are shown. In Figures 5.1 and 5.2 the signal is received from a simulation and in Figures 5.3 and 5.4 the signal is received from real data. For both situations the raw and the range corrected data are shown. When comparing the two signals the conclusion is that the simulated signal is similar to the real signal. The simulated signal is not as noisy as the real signal, but this can be changed by changing the atmospheric noise.

In Figure 5.5 the peak detector has been tested. Two different simulations with different atmospheric noise have been used. In the left image the atmospheric noise is less than in the right image. Different thresholds were tested and the lowest threshold for detection was set. In the left image the threshold was set to 6% above the standard deviation of the noise and in the right image the threshold was set to 21% above the standard deviation of the noise. In Figure 5.6 an image of the cloud concentration used in both images in Figure 5.5 is shown.

A simulated sweep is shown in Figure 5.7. In Figures 5.8 and 5.9 a simulation has been run with direct sunlight into the receiver. The electrical noise is clearly seen in the single waveform presented to the right in both figures. In Figures 5.10 and 5.11 a simulation has been run when the terrain is reflecting sunlight into the receiver. In the single shots presented to the right in the figures show less electrical noise compare to Figures 5.8 and 5.9.

(42)

30 Example of Results

Figure 5.1: Raw data from one simulated single shot in an empty atmosphere.

(43)

31

Figure 5.3: Raw data from one real single shot in an empty atmosphere, in com-parison with the simulated data in Figure 5.1.

Figure 5.4: A real, range corrected, single shot in an empty atmosphere, in com-parison with the simulated data in Figure 5.2.

(44)

Figure 5.5: Edge detection using different atmospheric noise. In the left image the atmospheric noise is 1% and the threshold is 6% above the noise level. In the right image the atmospheric noise is 5% and the threshold is 21% above the standard deviation of the noise level. The edges are marked with black dots.

(45)

33

Figure 5.7: A simulation of lidar signals when hitting a simulated cloud. The atmospherical noise has a standard deviation of 1 %.

(46)

Figure 5.8: A simulation when the sunlight is scattered into the receiver. A single shot is shown, when the laser is hitting the cloud, in the right image.

Figure 5.9: Same simulation as i Figure 5.8. The right image is showing a single shot in the empty atmosphere.

(47)

35

Figure 5.10: A simulation when reflections from the terrain are scattered into the receiver. A single shot is shown, when the laser is hitting the cloud, in the right image. Compare to Figure 5.8

Figure 5.11: Same simulation as i Figure 5.10. The right image is showing a single shot in the empty atmosphere. Compare to Figure 5.9

(48)

(49)

Chapter 6

Concluding Remarks

In this chapter a summary of the work will be presented and some conclusions and future directions will be discussed.

6.1 Summary

The goal of this thesis was to develop a simulation tool, present the data and test some detection algorithms on the received data.

During the work a simulation tool has been developed. A simple simulation script was already written and this was used as a base when a more complex simulation tool was developed. To make it easier to run simulations a graphical user interface was implemented. By using the GUI a lot of different parameters can be changed and the simulation can be executed without special knowledge of the program.

The next task was to present the results in a graphical user interface. After a simulation has been run the result can be shown in another GUI. In this GUI the detection algorithms can be applied and the result is shown as well.

Detection algorithms using single shots and the whole image have been imple-mented and tested on the simulated data.

6.2 Conclusions

The simulation tool developed is working well. When using the GUI a lot of different parameters can be changed to investigate the environmental and system influence on bio aerosol cloud detection. The concentration of the cloud discharged can be examined and the lidar signal as well. With this tool the smallest amount of particles to be detected can be examined. The detection algorithms work well on the simulated data. Depending on the atmospherical noise the detection limit is different. When the atmospherical noise is high the threshold must be set higher compared to a lower atmospheric noise. If the threshold is too low too much noise will be detected instead of the cloud. Subtraction of the background is one thing

(50)

38 Concluding Remarks

that needs to be done. In reality the laser will hit objects such as the ground, vegetation and buildings which needs to be subtracted from the signal. For real data the cloud might be divided into two or more different clouds which also need to be dealt with.

6.3 Future directions

For future projects it would be interesting to follow up some topics accordingly to • Take the laser hitting object such as the ground, vegetation and buildings

into consideration.

• Test the algorithms on real data.

• Implement the algorithms in C/C++ for a more rapid execution. With Matlab the simulation is quite slow and the presentation of the sweeps is slow as well.

• Simulation of two or more clouds. For a more realistic discharge the sim-ulation should be able to handle to simulate and also track two or more clouds.

• Calculate a warning area depending on the weather conditions.

• Implement different types of bio clouds and potential false clouds (dust, fog oil, diesel smoke etc.)

(51)

Bibliography

[1] V A Kovalev and W E Eichinger. Elastic lidar. Wiley interscience, USA, 1st edition, 2004. ISBN 0-471-20171-5.

[2] F Kullander, P Jonsson, P Wästerby, T Tjärnhage, G Olofsson, M Lindgren, and O Gustafsson. Ultraviolett laser- och detektionsteknik för kemisk och biologisk övervakning -förstudie. FOI, Oct 2007. FOI-R--2332--SE.

[3] L Persson. Invere Atmospheric Dispersion Modeling with Dynamical Bayesian Monte Carlo Methods. FOI CBRN Defence and Security, May 2008.

[4] A Pytel and J Kiusalaas. Engineering Mechanics Statics. Thomson Learning, Italy, 2nd edition, 2001. ISBN 1-86152-619-9.

[5] R Radke, S Andra, O Al-Kofahi, and B Roysam. Image change detection algorithms: A systematic survey. Aug 2004.

[6] O Steinvall. Laser system range calculations and the Lambert W function. Applied Optics Vol. 48, Issue 4, pp. B1-B7, Aug 2008.

[7] O Steinvall, R Persson, and O Gustafsson. Elastic Lidar Systems for Bio Aerosol Detection. FOI, Oct 2008. FOI-R--2455--SE.

[8] C Weitkamp and H Walther. Lidar: Range-Resolved Optical Remote Sensing of the Atmosphere. 2005. ISBN 0-387-40075-3.

[9] S R Yelisetty. Image change detection using wire-less sensor networks. India, Oct 2004. URL http://soar.wichita.edu/dspace/bitstream/10057/1185/1/t07061.pdf, Ac-cessed 17 November 2008.

(52)

(53)

Appendix A

Users Manual: Simulation

GUI

To make it easier to use the simulation a graphical user interface (GUI) has been developed. The GUI is a link between the code and the user. In this Appendix some screenshots from the simulation will be presented and instruction of how to use the program is found here.

A.1 Cloud Concentration

In the GUI the user can change a lot of input parameters for the simulation. The first step is to simulate the cloud concentration with the Umeå model. Parameters to change can be seen in Figure A.1. The A and B parameters is detemine the wind strength in the start moment of the simulation in the x- and y-directions. When all parameters has been chosen the cloud concentration, Ncloud, is simulated

by clicking on the button Get Cloud. The simulated cloud concentration can be visualised in a colour plot by using the button Show Cloud Conc where the concentration is given in particles per litre.

A.2 The Lidar Signal

In the next step the lidar signals from scanning the atmosphere with the simu-lated cloud are simusimu-lated. A lot of parameters have to be determined before the simulation can be made. A screen shot from the GUI can be seen in Figure A.2.

Laser

The first group of parameters is consisting of laser parameters. Here all the laser parameters can be changed, if not, the default value will be used. The parameters to change are the pulse energy and the pulse length.

(54)

42 Users Manual: Simulation GUI

Figure A.1: A part of the GUI where the cloud concentration is simulated.

Receiver

The receiver parameters are influence on the noise calculations. Here the diameter and the focal lenght of the receiver can be chosen. The overlap range and the range where the overlap has reached 50 % is also changed here.

Optical Background

The background noise is consisting of reflected sun light and terrain reflectivity into the receiver.

Receiver Noise

In those boxes the calculated NEP will be presented after the simulation has run.

Atmosphere

Here the visibility range can be changed. In the yellow boxes the calculated backscatter coefficients will be shown after the simulation has been run.

B Cloud

If choosing Cloudsim the received signal is pointing at the simulated cloud from Umeå. At the moment the particle used is BG. By changing the the popup menu to Optional coefficients the coefficients βc and αc can be set by the user. The

(55)

A.2 The Lidar Signal 43

Figure A.2: The simulation GUI.

other choice, Empty Atmosphere, an empty atmosphere will be simulated without any cloud in the field of view.

Source

The coordinates where the discharge starts can be set here and also the radius and standard deviation of the particles.

Search Parameters

Here the angel resolution and the range can be set. At the moment the scanning angle is set to 30 and cannot be changed.

The Saved Matrixes

The simulated data are stored in mat-files and can be received by writing load in the Matlab command window and thereafter the wanted matrix. The raw data matrix is named P_raw.mat and the range corrected data matrix is named P_corrected.mat.

(56)

Appendix B

Users Manual: Presentation

GUI

In this chapter the presentation graphical user interface will be presented. The presentation GUI can be seen in Figure B.1.

(57)

B.1 Explanation of the Action of the Buttons 45

B.1 Explanation of the Action of the Buttons

In this section all buttons in the presentation GUI will be explained.

Show Sweeps

A beforehand simulated cloud is presented in the upper left square. All sweeps are showing after each other, the pulse repetition frequency can be set by using the box next to the Show Sweeps button. Below the Show Sweeps button the time for each sweep is presented. By inserting a number in the box above the colourbar next to the upper left square the colourmap span is changed the next time the

Show Sweeps is used. The x- and y-axes are showing the range in meters and

the colour is indicating the amount of backscattered light. Red is indicating high backscattering and blue is indicating low backscattering.

Mouse Click

When clicking on the Mouse Click button the mouse cursor is transformed to a cross, see Figure B.2. By clicking in the PPI-image in the upper left square the chosen single waveform will be shown in the upper right square. The waveform number is shown above the square, the x-axis is showing the range in meters and the y-axis is showing the signal amplitude. By using the Mouse Click button again another waveform can be chosen. To the right of the Mouse Click button there are two radiobuttons. These buttons will chose which data to show, the range corrected or the raw data.

Figure B.2: The mouse cursor is changing when clicking on the Mouse Click button and a waveform can be chosen by clicking in the PPI-image.

Waterfall

The Waterfall button is making a waterfall plot of the present sweep showed in the upper left square. The waterfall plot will be shown in the lower right square,

(58)

46 Users Manual: Presentation GUI

see Figure B.3.

Figure B.3: A waterfall plot is shown when activation of the Waterfall button.

Show Waveform

Instead of using the mouse to chose a single waveform to be shown the waveform number can be typed into the box next to the Show Waveform button. The chosen waveform number will then be shown in the upper right square when the Show

Waveform button is activated.

Peak Find

This button will find the right and left edge of the cloud and mark it with white stars in an PPI-image in the lower left square, see Figure B.4. The threshold used is a certain percentage above the noise level. By changing the number in the box next to the button the threshold value is changed.

Centre of Mass

When using this button the centre of mass of the cloud is calculated and marked with a green star in a PPI-image in the lower left square. The coordinates for the centre of mass is printed in the boxes next to the button.

Track Cloud

If this button is active the cloud in the scene will be tracked by a narrower sweep, see Figure B.5 The threshold used can be chosen by changing the threshold box.

(59)

Figure B.4: Showing the right and left edge found by the peak detection algorithm.

Sweep + Conc

If wanted, both the lidar signals and the concentrations can be illustrated at the same time. By using the button Sweep + Conc the lidar signals will be shown in the upper left square and the concentration plot will be shown in the lower left square, the concentration is given in particles per litre as can be seen in Figure B.7. The colourmap for the concentration image can be changed by inserting a value in the box above the colourbar and pressing enter.

Stop

With this button the cloud tracking and Show Sweeps can be stopped.

Clear All

(60)

48 Users Manual: Presentation GUI

Figure B.5: When the botton Track Cloud is chosen the sweeps will be narrower.

Figure B.6: Changing the colourmap span determines the cloud visibility in the image.

(61)

Figure B.7: The upper image is showing the lidar signal and the lower image is showing the concentration in particles per litre.

A Simulation Model for Detection and Tracking Bio Aerosol Clouds using Elastic Lidar

Institutionen för systemteknik

Department of Electrical Engineering

Examensarbete

A Simulation Model for Detection and Tracking Bio

Aerosol Clouds using Elastic Lidar

A Simulation Model for Detection and Tracking Bio

Aerosol Clouds using Elastic Lidar

Examensarbete utfört i Reglerteknik

vid Tekniska högskolan i Linköping

av

Abstract

Sammanfattning

Acknowledgments

Contents

Chapter 1

Introduction

1.1

Problem Description

1.2

Purpose

1.3

Outline

Chapter 2

Lidar Theory

2.1

The Lidar Equation

2.1.1

Power of the Laser P

2.1.2

Optical Efficiency η

2.1.3

Range Resolution

2.1.4

The Overlap Function O(R)

2.1.5

Backscatter Coefficient β

(R, λ)

2.1.6

Atmospheric Transmisson T

(R, λ)

2.1.7

Noise n(R)

Chapter 3

Simulation of Bio Aerosol

Discharges

3.1

The Simulation Model

3.1.1

The Overlap Function

3.1.2

Backscatter

3.1.3

Transmission

3.1.4

Noise

3.2

Calculation of the Range in Clear Atmosphere

3.3

Cloud Concentration

Chapter 4

Detection and Tracking of

Simulated Bio Aerosol

clouds

4.1

Pre Processing

4.1.1

Range Correction

4.1.2

Moving Average

4.2

Change Detection

4.3

Single Waveform Analysis

4.3.1

Finding Maxima

4.3.2

Simple Peak Detector

4.3.3

Removing False Peaks

_c

_c