sensors Article Smelling Nano Aerial Vehicle for Gas Source Localization and Mapping Javier Burgués 1,2, * , Victor Hernández 3 , Achim J. Lilienthal 3 and Santiago Marco 1,2 1 Institute for Bioengineering of Catalonia (IBEC), The Barcelona Institute of Science and Technology, Baldiri Reixac 10-12, 08028 Barcelona, Spain; smarco@ibecbarcelona.eu 2 Department of Electronics and Biomedical Engineering, Universitat de Barcelona, Marti i Franqués 1, 08028 Barcelona, Spain 3 AASS Mobile Robot Olfaction Lab, Örebro University, SE 70182 Örebro, Sweden; victor.hernandez@oru.se (V.H.); achim@lilienthals.de (A.J.L.) * Correspondence: jburgues@ibecbarcelona.eu; Tel.: +34-934-029-070  Received: 11 December 2018; Accepted: 22 January 2019; Published: 24 January 2019  Abstract: This paper describes the development and validation of the currently smallest aerial platform with olfaction capabilities. The developed Smelling Nano Aerial Vehicle (SNAV) is based on a lightweight commercial nano-quadcopter (27 g) equipped with a custom gas sensing board that can host up to two in situ metal oxide semiconductor (MOX) gas sensors. Due to its small form-factor, the SNAV is not a hazard for humans, enabling its use in public areas or inside buildings. It can autonomously carry out gas sensing missions of hazardous environments inaccessible to terrestrial robots and bigger drones, for example searching for victims and hazardous gas leaks inside pockets that form within the wreckage of collapsed buildings in the aftermath of an earthquake or explosion. The first contribution of this work is assessing the impact of the nano-propellers on the MOX sensor signals at different distances to a gas source. A second contribution is adapting the ‘bout’ detection algorithm, proposed by Schmuker et al. (2016) to extract specific features from the derivative of the MOX sensor response, for real-time operation. The third and main contribution is the experimental validation of the SNAV for gas source localization (GSL) and mapping in a large indoor environment (160 m2 ) with a gas source placed in challenging positions for the drone, for example hidden in the ceiling of the room or inside a power outlet box. Two GSL strategies are compared, one based on the instantaneous gas sensor response and the other one based on the bout frequency. From the measurements collected (in motion) along a predefined sweeping path we built (in less than 3 min) a 3D map of the gas distribution and identified the most likely source location. Using the bout frequency yielded on average a higher localization accuracy than using the instantaneous gas sensor response (1.38 m versus 2.05 m error), however accurate tuning of an additional parameter (the noise threshold) is required in the former case. The main conclusion of this paper is that a nano-drone has the potential to perform gas sensing tasks in complex environments. Keywords: robotics; signal processing; electronics; gas source localization; gas distribution mapping; gas sensors; drone; UAV; MOX sensor; quadcopter 1. Introduction Thanks to recent advances in micro-technology, manufacturers of unmanned aerial vehicles (UAVs), or drones, have been able to develop miniaturized flying platforms; with insect-sized aircrafts expected in the future [1] (Figure 1). A micro-UAV (MAV or µUAV) has a length between 15 cm and 100 cm and a weight between 50 g and 2 kg [2]. A nano air vehicle (NAV) or nano-drone is extremely small, with a wingspan lower than 15 cm, and weighs less than 50 g [2]. If compared to piloted aircrafts Sensors 2019, 19, 478; doi:10.3390/s19030478 www.mdpi.com/journal/sensors Sensors 2019, 19, 478 2 of 25 or larger Sensors UAVs, 2018, 18, x FOR PEER REVIEW MAVs and NAVs can fly at low altitudes (i.e., below 150–200 m) over small 2 of 26 geographic or site-specific areas on a real-time basis at affordable operational costs [3]. MAVs piloted aircrafts or larger UAVs, MAVs and NAVs can fly at low altitudes (i.e., below 150–200 m) equipped with gas detection systems over small and/or geographic or sampling site-specific bags areas have been already on a real-time basis at used in the affordable fields ofcosts operational environmental [3]. monitoringMAVs[4–9], volcanic equipped withgas gas sampling [10–14], detection systems localization and/or sampling of fugitive bags emissions have been [15,16], already used in theearly fire detectionfields of environmental [17,18], monitoring[19–21], precision agriculture [4–9], volcanic landfillgas sampling [10–14], monitoring [22–24],localization of fugitive[25] and disaster response emissions[15,16], early fire detection [17,18], precision agriculture [19–21], landfill monitoring [22– mine blasting [26], among others [27,28]. 24], disaster response [25] and mine blasting [26], among others [27,28]. Figure 1.Figure 1. Overview Overview of theofUAV the UAV landscape,from landscape, frominsect-sized insect-sized drones to military drones aircrafts, to military classified aircrafts, classified according to the approximate weight and size. The graphic shows the large range of UAV sizes, which according to the approximate weight and size. The graphic shows the large range of UAV sizes, spans seven orders of magnitude. which spans seven orders of magnitude. The tiny form-factor and maneuverability of NAVs allow sensing of hazardous environments The inaccessible tiny form-factor androbots to terrestrial maneuverability of NAVs and bigger drones, can flyallow sensing over areas being of hazardous unobserved and environments are not inaccessible to terrestrial a hazard for humans, robots andtheir enabling bigger use indrones, can fly public areas over buildings. or inside areas being unobserved Providing and are not a NAV with olfaction capabilities is now possible due to miniaturization and low-cost fabrication a hazard for humans, enabling their use in public areas or inside buildings. Providing a NAV with of gas sensors. olfactionAn odor-sensitive capabilities nano-drone is now possiblecan be dueused to in a myriad of applications miniaturization that rangefabrication and low-cost from environmental of gas sensors. monitoring to search and rescue, leak detection, chemical, biological, radiological and nuclear An odor-sensitive nano-drone can be used in a myriad of applications that range from (CBRN) defense, explosive finding, among others. For example, in the aftermath of an earthquake or environmental monitoring to search explosion and rescue, it is important leak for to search detection, victims andchemical, hazardous biological, radiological gas leaks inside and pockets that nuclear form within (CBRN) defense, explosive the wreckage finding, among of collapsed others. buildings. A For example, nano-drone innavigate could the aftermath of an earthquake such scenarios or explosion much faster than a terrestrial robot, passing through confined spaces that preclude human entry, it is important to search for victims and hazardous gas leaks inside pockets that form within the evading obstacles or wreckagelarge gaps and sampling the space in three dimensions (3D). of collapsed buildings. A nano-drone could navigate such scenarios much faster than a terrestrial robot, passing through confined 1.1. Related Work on Gas-Sensitive Nanodronesspaces that preclude human entry, evading obstacles or large gaps and sampling the space in three dimensions (3D). Two experimental works [29,30] already explored the viability of nano-drones for gas sensing tasks. Rossi et al. [29] performed preliminary indoor experiments using a CrazyFlie 2.0 nano-drone 1.1. Related Work on Gas-Sensitive Nanodrones equipped with a metal-oxide semiconductor (MOX) gas sensor. The authors found that the air drawn Twoaround the airframe experimental worksstrongly affected [29,30] alreadythe explored sensor response, resulting the viability of innano-drones useless signals. forThey gas sensing evaluated several mechanical solutions to keep the sensor out of the region of airflow created by the tasks. Rossi et al. [29] performed preliminary indoor experiments using a CrazyFlie 2.0 nano-drone nano-drone, but the drone became uncontrollable because of inertia problems. The adopted solution equippedwas with a metal-oxide to operate the dronesemiconductor (MOX) gas in the so-called “butterfly” sensor. mode, Theaauthors in which found human pilot that lands the the air drawn drone around the airframe strongly affected the sensor response, resulting in useless signals. in the proximity of the source and halts the motors to take a measurement. In this way, the sensor They evaluated signals are not several mechanical affected to solutions by keep the rotors the but, at the sensor outsame time, of the the 3Dofsensing region airflow capabilities created of bythe thedrone nano-drone, are not but the drone used, and became the approach might uncontrollable not scale because well to problems. of inertia large scenarios. The adopted solution was to operate Fahad et al. [30] equipped the same nano-drone with a chemically sensitive field effect transistor the drone in the so-called “butterfly” mode, in which a human pilot lands the drone in the proximity of (CS-FET) sensor for hydrogen (H2) detection. The test environment was a chemical hood in which H2 the source and halts the motors to take a measurement. In this way, the sensor signals are not affected by the rotors but, at the same time, the 3D sensing capabilities of the drone are not used, and the approach might not scale well to large scenarios. Fahad et al. [30] equipped the same nano-drone with a chemically sensitive field effect transistor (CS-FET) sensor for hydrogen (H2 ) detection. The test environment was a chemical hood in which H2 was released from the upper part of the hood and the drone ascended from the bottom of the hood to Sensors 2019, 19, 478 3 of 25 the area near the source (h = 60 cm), aided by high tension strings. The sensor response increased as the drone approached the source, reaching its maximum value after hovering (i.e., levitating) near the source for 40 s. The drone was merely used as a proof-of-concept demonstration of the proprietary gas sensor developed by the authors. The above works suggest that a nano-drone might be used for gas source localization (GSL), however the experimental scenarios were extremely simple. 1.2. Experimental Evaluation of Gas-Sensitive Nanodrones The main problems for performing large-scale experiments in complex environments using nano-drones are related to the limited on-board resources and difficulty to control the platform due to inertia and stability issues. Taking as an example the CrazyFlie 2.0 (CF2) quadcopter, the tiny 240 mAh battery delivers power for up to 7 min of flight and 15 grams of payload, which means that only lightweight and power-efficient sensors can be mounted on board. Self-localization and obstacle avoidance—required for autonomous navigation—are hard to accomplish because laser scanners, for example, are too heavy for the nano-drone payload. Autonomous hovering of a CF2 equipped only with a tiny RGB camera has been achieved in indoor experiments [31], although the camera consumes all available payload and reduces the flight time to 3.5 min. The Global Positioning System (GPS) can be used for localization outdoors where, however, nano-drones can often not be controlled stably due to their low inertia and strong wind. Navigation in indoor areas can be achieved through motion capture systems (MOCAPs) [32] or radio frequency (RF)-based systems [33]. MOCAPs offer high accuracy (1 mm error) but are expensive, typically only cover small volumes and require line-of-sight (LOS) between the cameras and the drone. RF-based systems are cheaper, have a larger coverage area, do not necessarily require LOS but are less accurate (10 cm error). In many realistic scenarios, deploying an external localization system might not be possible (e.g., in a disaster situation) and the drone would have to navigate autonomously or remotely controlled by an operator. Gas sensing tasks are also subject to additional constraints, as they must be executed in the short time limited by the battery capacity and relying exclusively on one or two chemical sensors. It should be noted that most research on GSL is based on terrestrial robots, which can be running for hours, perform long measurements of 1–2 min at each sampling location and possibly use selective sensors (e.g., TDLAS, OGI cameras, e-noses) and anemometers. Large drones can fly for 20–25 min and be equipped with the same technology as terrestrial robots. Nano-drones are therefore subject to unprecedented constraints because a stop-sense-go strategy would only allow for at most 10 measurements (of 30 s each), the limited number of sensors hinder rejecting chemical interferences and the absence of anemometry prevents assessing the wind direction, which is a key parameter for GSL. During its operation, the drone can also not fly too fast or the relevant structures of the chemical plume may become blurred due to the slow response time of the sensors [34,35]. Lilienthal et al. [35] pointed out that the gas distribution mapped by a terrestrial robot may be slightly shifted as compared to the real distribution, due to the memory effect of MOX sensors. 1.3. Gas Source Localization Gas source localization (GSL) is a key task for gas-sensitive robots that consists in identifying the point of release of a hazardous gas. GSL strategies can be divided into three groups: reactive plume tracking (bioinspired), plume modelling and gas distribution mapping (GDM) strategies [36] (Figure 2). Bioinspired algorithms attempt to track the gas plume along its entire length, mimicking the excellent odor plume tracking capabilities of some flying insects. At this time it is unclear whether bioinspired reactive behaviours have better performance than other approaches based on statistical inference from cumulative readings [37–39]. According to [11], the bioinspired reactive behaviors that researchers have implemented on mobile robots are modelled too simple to cope with complex environments, gas sensors used are too slow to resolve plume features in a milliseconds scale and mobile robots are not agile enough for performing insect-like reactive movements. Besides of that, bioinspired algorithms often require wind measurements (anemometry) and real-time obstacle detection, which hinders their use in nano-drones. Sensors 2019, 19, 478 4 of 25 Plume modelling algorithms [40–43] assume a mathematical model for the plume, such as Gaussian shaped Sensors 2018, 18,plumes x FOR PEER [44] or filament/particle based models [45–47], and use local measurements REVIEW 4 of 26 of concentration and wind to fit the model and estimate the source location, which is usually a parameter of the model.Plume The modelling practicalalgorithms [40–43] applicability assume modelling of plume a mathematical modelisfor methods the plume, limited because suchtheyas tend to makeGaussian overly shaped plumes [44] simplifying or filament/particle assumptions based (e.g., that themodels wind[45–47], field isand use local stable, measurements spatially uniform and of concentration and wind to fit the model and estimate the source location, which is usually a measurable), often require a-priori information such as the source release rate in Gaussian models [40], parameter of the model. The practical applicability of plume modelling methods is limited because or are they sensitive to meta-parameters such as the odor detection threshold in filament-based models [41] tend to make overly simplifying assumptions (e.g. that the wind field is stable, spatially uniform or theand probability measurable), of particle encounter often require as ainformation a-priori function ofsuch distance as thetosource the source releaseinrate particle models [42]. in Gaussian The Gaussian models [40], model or arealso assumes sensitive that the exploration to meta-parameters such as the area odordoes not contain detection thresholdobstacles in filament- or walls whichbased couldmodels otherwise [41] ordistort the plume. the probability Further, of particle long time-averaging encounter as a function ofmightdistancebetorequired the source to inobserve particle models [42]. The Gaussian model also assumes that the exploration a Gaussian plume [48,49] or to estimate the particle density in a certain region, which slows down area does not contain plumeobstacles modelling or walls which could approaches otherwise distort the plume. Further, long time-averaging might be for GSL. required to observe a Gaussian plume [48,49] or to estimate the particle density in a certain region, GDM approaches use sensor measurements to first build a map of the gas distribution in the which slows down plume modelling approaches for GSL. environment, which is then used to estimate the source location. Maps reflecting the instantaneous GDM approaches use sensor measurements to first build a map of the gas distribution in the concentration environment, [50],which the mean is thenconcentration used to estimate[35,38], the location. the source varianceMapsof the concentration reflecting [38,51] or the the instantaneous number of odor hits concentration [50],(which the mean areconcentration over-threshold segments [35,38], in the the variance sensor of the response) concentration [52]orhave [38,51] the been successfully numberused of odorfor GSL. To build hits (which are aover-threshold gas distribution map, the segments path in the of the sensor robot should response) [52] haveroughly been cover successfully the entire search usedarea,for GSL. To moving typically build a gas distribution along map, thetrajectory a predefined path of the[35,38,50], robot should roughly cover although adaptative the entire search area, typically moving along a predefined trajectory approaches have been proposed [53]. When the map is based on statistical properties of the gas[35,38,50], although adaptative approaches have been proposed [53]. When the map is based on statistical properties of the gas distribution (e.g., mean or variance), long measurements (typically 30 s and more) are often carried distribution (e.g., mean or variance), long measurements (typically 30 s and more) are often carried out at each sampling location. This stop-and-sense strategy is not suitable for UAVs, particularly for out at each sampling location. This stop-and-sense strategy is not suitable for UAVs, particularly for nano-drones, sincesince nano-drones, any hovering any hovering stop stop quickly drains quickly the battery drains [15].[15]. the battery GDM GDM algorithms algorithms areare less efficient less than bioinspired efficient than bioinspired and plume modelling algorithms (in terms of distance travelled by the but and plume modelling algorithms (in terms of distance travelled by the robot) do notrobot) rely on butunrealistic do not rely assumptions nor require on unrealistic assumptionswindnor information require wind or a-priori parameters. information or a-priori They are parameters. compatible with They slower are robots compatible andwith slower sensors androbots and sensorsgas the resulting anddistribution the resultingmap gas distribution can be used for other map can bebeyond purposes used for GSL.other purposes beyond GSL. 2. Gas FigureFigure 2. source localization Gas source strategies. localization strategies. (left) Reactiveplume (left) Reactive plume tracking; tracking; (center) (center) Plume Plume modelling; modelling; (right)(right) Map-based. Map-based. MOXMOX sensors are are sensors probably probably the most the mostsuitable gassensing suitable gas sensingtechnology technology for nano-drones for nano-drones due toduethe to the reduced size (few reduced mm2mm size (few ), low power 2), low powerconsumption consumptionofofsome somemodels models (few mW)[54] (few mW) [54]andandsimplicity simplicity of of the the conditioning electronics. Lilienthal et al. [55] observed that the maximum conditioning electronics. Lilienthal et al. [55] observed that the maximum response of a MOX sensor response of a MOX sensor often corresponds often corresponds to the approximate to the approximate location of location a gas of a gasifsource source if the readings the sensor sensor readings are are acquired in acquired in motion. Such a correlation was never observed if the concentration measurements were motion. Such a correlation was never observed if the concentration measurements were collected collected with a stop-sense-go strategy. This behaviour, which was previously reported by Atema et with a stop-sense-go strategy. This behaviour, which was previously reported by Atema et al. [56] al. [56] and confirmed by Farah and Duckett [57], was attributed to the long recovery time of MOX and confirmed sensors. If by Farah a MOX and Duckett sensor is exposed [57], wasconsecutive to two attributed gas to the long recovery patches, the response timeto of theMOXsecondsensors. If a MOX sensor stimulus willis occur exposedwhen to the twosensor consecutive has not gas yet patches, recoveredthe from response the firsttoexposure. the second Thestimulus overall will occur response when thetosensor has not the second patch yetwill recovered be higherfrom than the first if the exposure. sensor had been The overall already response fully recovered to from the second patch the willfirst be patch. higherSincethantheif local densityhad the sensor of gas patches been tendsfully already to berecovered higher nearfromthe source, thepatch. the first encounterSince the rate with gas patches is higher if the robot is moving. Thus, it is plausible local density of gas patches tends to be higher near the source, the encounter rate with gas patches that Rossi et al. [36] and Luoifetthe is higher al. robot [50] obtained is moving. a good Thus, estimate of source it is plausible location that Rossi in et outdoor experiments al. [36] and Luo et al.using [50] the obtained instantaneous response of a MOX sensor mounted on a micro drone ( 800 g). In the latter case, the a good estimate of source location in outdoor experiments using the instantaneous response of a MOX sensor mounted on a micro drone (≈800 g). In the latter case, the authors were able to build a 3D gas Sensors 2019, 19, 478 5 of 25 distribution map of a relatively large outdoor environment (10 × 16 m2 ) and localize the source in less than 10 min. For that, they used a predefined flight path consisting of two 2D rectangular sweepings at different heights (0.3 m and 1 m), without stopping at predefined locations for measuring. The cell with maximum value of the gas distribution map coincided approximately with the true source location. The results of Luo et al. [50] also show that the density of odor hits in a map is correlated to the source location. An odor hit is typically declared when the instantaneous concentration exceeds a certain threshold. Although the term ‘odor hit’ is widely used in the literature to indicate contact between the gas and the sensor, strictly speaking it would be more precise to refer to these events as ‘gas hits’ or ‘plume hits’ since most gases being released in reported experiments are odorless for humans. Nonetheless, Thomas Lochmatter [11] argues that the definition of odor as a gas that humans can smell can be extended to robots. In this sense, terms such as “odor” can also be used instead of “gas”. Odor hits are supposedly caused by contact with individual patches of the plume and there are indications that insects use similar features to orient rapidly in turbulent plumes [58,59]. Well-known GSL algorithms such as Pang and Farrell’s method [41] or Infotaxis [42] model the plume as a sequence of chemical filaments/particles and use odor hits to localize the source. Detecting odor hits with a MOX sensor is a challenging task due to the long recovery time (in the order of 10–30 s) compared to the temporal resolution of the chemical stimuli (in the order of ms). Several research groups attempted to improve the response time of MOX sensors by novel hardware designs [60,61] or by signal processing, using inverse dynamical models [34,62,63], artificial neural networks [64] or extracting specific features [65]. Hardware methods are not appropriate for nano-drones, as they usually increase the size, weight and power consumption of the system. Signal processing methods are more suitable. Schmuker et al. [65] proposed a method to extract short time-scale features (called ‘bouts’) from the derivative of the MOX sensor response that could be caused by contact with individual filaments of the plume. In wind tunnel experiments, the authors found that the frequency of these ‘bouts’ (as detected with MOX sensors) is strongly correlated to the distance of a gas source: the higher the bout frequency, the closer the sensor to the gas source. The proposed algorithm uses a threshold to filter out low-amplitude bouts—produced by sensor noise—that would otherwise lead to meaningless correlations. It was also found that the variance of the bout frequency (measured across multiple trials) indicates whether the detector is in the plume centerline (low variance) or slightly lateral from it (high variance). These features suggest a plume-tracking GSL strategy in which the robot first tries to locate the plume centerline by monitoring the variance of the bout frequency and then approaches the source by moving in the direction of increasing bout frequency. The advantages of this method are that anemometry is not required and, since bouts are detected in the derivative of the signal, the algorithm is not very sensitive to changes in the background concentration or to differences between individual gas sensors. The sensitivity of the algorithm to the threshold used to discard noise-induced bouts has not been studied yet, but it might have a large impact on the results. Since this method has not been experimentally validated beyond a wind tunnel, it remains to be shown if meaningful gradients of both bout frequency and its variance can be obtained in real scenarios. 1.4. Proposed Smelling Nano Aerial Vehicle (SNAV) In this work, we propose a Crazyflie 2.0 nano-drone equipped with a MOX sensor for gas source localization in large indoor environments. We calibrate the sensor to compensate the non-linear response, obtain measurements in concentration units and to estimate the limit of detection (LOD). We assess the impact of the propellers on the MOX sensor signals at different distances of a chemical source. We then compare two GSL strategies, one based on the instantaneous response and the other one based on the bout frequency in two experiments where the source is placed challenging positions for the drone. We show that proper selection of the bout amplitude threshold is critical for good localization performance. We also demonstrate that a 3D gas distribution map of an environment of 160 m2 can be built in less than 3 min using the proposed platform and the source can be accurately localized from the map. Sensors 2019, 19, 478 6 of 25 2.Sensors 2018, 18, Materials andx FOR PEER REVIEW Methods 6 of 26 2. Materials 2.1. Nano-Droneand Methods and Gas Sensors 2.1.Among Nano-Drone all commercial and Gas Sensors nano-drones, we selected the CF2 (Bitcraze AB, Malmö, Sweden) due to its low cost, reduced dimensions (10 × 10 cm) and open hardware/software architecture. Weighing Among all commercial nano-drones, we selected the CF2 (Bitcraze AB, Malmö, Sweden) due to only 27 g, it has a maximum recommended payload of 15 g and is capable of up to 7 min of continuous its low cost, reduced dimensions (10 × 10 cm) and open hardware/software architecture. Weighing flight. The main microcontroller (µC) is an ARM 32-bit STM32F405 Cortex-M4, which runs an open only 27 g, it has a maximum recommended payload of 15 g and is capable of up to 7 minutes of source real-time operating system (FreeRTOS). The CF2 communicates with a ground station (PC with continuous flight. The main microcontroller (μC) is an ARM 32-bit STM32F405 Cortex-M4, which USBrunsradio antenna) an open source over the 2.4operating real-time GHz ISMsystem radio band in up toThe (FreeRTOS). 1 km CF2range line-of-sight communicates (LOS). with The CF2 a ground periodically sends over this link internal variables of the system station (PC with USB radio antenna) over the 2.4 GHz ISM radio band in up to 1 km range line-of- (sensor measurements, position, rotor sightspeed, (LOS).battery The CF2 level, user-defined periodically variables, sends over this etc.) andinternal link receives commands variables of thefrom the base system station (sensor such as waypoints measurements, or parameter position, updates. rotor speed, An expansion battery port—accessible level, user-defined variables,frometc.)topand andreceives bottom of the drone—provides commands from the the baseuser with station access such to certainorµC as waypoints pins (4 xupdates. parameter GPIO, 3An x PWM), expansion power port— lines (GND, VCC 3.0 V, VBAT 3.7 V) and communication buses (I2C, accessible from top and bottom of the drone—provides the user with access to certain μC pins (4 x SPI, 2 x UART). New hardware compatible with 3.0V GPIO, 3 x PWM), powerlogiclines can (GND, be easily VCCintegrated 3.0 V, VBAT into 3.7 theV)platform by soldering buses and communication it into(I2C, an expansion SPI, 2 x UART). board, calledNew hardware a deck, which compatible with 3.0V can be connected logic to the can be easily expansion port. integrated An installed into theisplatform deck by automatically soldering it into an expansion board, called a deck, which can be connected detected and initialized by the CF2 at startup, without having to modify the stock firmware. Only the to the expansion port. An installed deck driverdeckneeds is automatically to be programmed, detected and initialized a separate piece ofbycode the that CF2 specifies at startup, thewithout having of functionality to the modify deck the stock firmware. in FreeRTOS language. Only the deck driver needs to be programmed, a separate piece of code thatAspecifies the functionality custom deck (i.e., a printed of the deck in circuit FreeRTOS board), named language. the MOX deck, was developed to interface A custom deck (i.e., a printed circuit two MOX gas sensors to the CF2 (Figure 3). The deck contains board), named the MOXtwo deck, was developed sockets for 4-pinto interface Taguchi-type two MOX gas sensors to the CF2 (Figure 3). The deck contains (TGS) gas sensors, a temperature/humidity sensor (SHT25, Sensirion AG, Stäfa, Switzerland), two sockets for 4-pin Taguchi-type (TGS) gas sensors, a temperature/humidity sensor (SHT25, Sensirion AG, Stäfa, Switzerland), a dual- a dual-channel digital potentiometer (AD5242BRUZ1M, Analog Devices, Norwood, MA, USA) and channel digital potentiometer (AD5242BRUZ1M, Analog Devices, Norwood, MA, USA) and two two MOSFET p-type transistors (NX2301P, NEXPERIA, Nimega, NL). We selected the TGS 8100 sensor MOSFET p-type transistors (NX2301P, NEXPERIA, Nimega, NL). We selected the TGS 8100 sensor (Figaro Engineering Inc., Osaka, Japan) due to its compatibility with 3.0 V logic, power consumption of (Figaro Engineering Inc., Osaka, Japan) due to its compatibility with 3.0 V logic, power consumption only 15 mW (the lowest in the market as of June 2016) and miniaturized form factor (MEMS). Since the of only 15 mW (the lowest in the market as of June 2016) and miniaturized form factor (MEMS). Since sensor heater the sensor usesuses heater 1.8V,1.8V, we included we included twotwotransistors transistors (one (oneper persensor) sensor)totoreduce reducethe theapplied power by applied power means of pulse by means width of pulse width modulation modulation (PWM). (PWM).The TheMOXMOX read-out circuit(Figure read-out circuit (Figure4)4) is is a voltage a voltage divider divider connected to the µC’s analog-to-digital converter (ADC). The voltage connected to the μC’s analog-to-digital converter (ADC). The voltage divider is powered at 3.0 V anddivider is powered at 3.0 V and thethe load resistor (RLL) can be set dynamically by the potentiometer (from 60 Ω to 1 MΩ in steps of 3.9 of load resistor (R ) can be set dynamically by the potentiometer (from 60 Ω to 1 MΩ in steps 3.9kΩ). kΩ).InInthethe current current work, work, wewe only only used used oneone of ofthethetwotwo sensors sensors and and RL Rwas L was fixed fixed to 70to kΩ. 70 kΩ.ThisThis loadload resistor resistor value is selected to operate the voltage divider near its mid-range, where the sensitivity is is value is selected to operate the voltage divider near its mid-range, where the sensitivity maximum maximum(according (accordingto tothe the expected concentrationsfor expected concentrations forthethegasgassource sourcedescribed described in in Section Section 2.2). 2.2). Figure3.3.The Figure TheCrazyFlie CrazyFlie 2.0 2.0 equipped withthe equipped with theMOX MOXdeckdeckand andthe the UWB UWB tagtag (center) (center) gets gets its 3D its 3D position position from fromananexternal externallocalization localization system composedof system composed ofsix sixultra-wide ultra-wideband band anchors anchors (left). (left). TheThe location location andand sensordata sensor dataare arecommunicated communicated to the ground ground station station(right) (right)over overthethe2.4 2.4GHz GHz ISMISMband. band. The Theinitialization initializationtask task of of the MOX MOXdeck deckdriver driverconfigures configures the the PWM, PWM, initializes initializes thethe SHT25 SHT25 sensor, sensor, setsthe sets thewiper wiperposition positionof of both both channels channels of ofthe thepotentiometer potentiometerandandadds thethe adds MOX MOX readout registers readout to to registers the list of variables that are continuously logged and transmitted to the base station. The goal of PWM Sensors 2019, 19, 478 7 of 25 Sensors the 2018, list of 18, x FORthat variables PEER REVIEW are 7 of 26 continuously logged and transmitted to the base station. The goal of PWM is to convert from 3.0 V to 1.8 V required by the MOX heater. For that, the PWM frequency is set to is toSensors convert from 3.0 V to 1.8 V required by the MOX heater. For that, the PWM frequency is set to 8.4 KHz and the 2018, 18, duty cycleREVIEW x FOR PEER (DC) to 36% according to 7 of 26 8.4 KHz and the duty cycle (DC) to 36% according to is to convert from 3.0 V to 1.8 V required by the2MOX heater.2For that, the PWM frequency is set to DC𝐷𝐶 = 𝑃 = Vavg𝑉/R/𝑅 Pavg = 1.8 1.8 = 36%, (1) 8.4 KHz and the duty cycle (DC) to 36% according Ppeak 2 to 2 36%, (1) 𝑃 𝑉 /R /𝑅 3.03.0 Vpeak 𝑃 𝑉 /𝑅 1.8 𝐷𝐶 36%, (1) where 𝑃 is is wherePavg the the averagepower average powerdelivered𝑃 totothe delivered the /𝑅 P3.0 𝑉 sensor, sensor, 𝑃 is is peak thethe peak peak power power ofof thePWM the PWMsignal signal 2 /R (Joule’s first law combined with Ohm’s law) is used to convert from power to voltage. and 𝑃 𝑉 /𝑅 (Joule’s first law combined with Ohm’s law) is used to convert from power to voltage. and P = V where 𝑃 is the average power delivered to the sensor, 𝑃 is the peak power of the PWM signal The Theduty dutycycle cycleisisthe thefraction fractionofoftime timethat thatthe thetransistor transistordelivers deliverspower power totothe sensor. the sensor. and 𝑃 𝑉 /𝑅 (Joule’s first law combined with Ohm’s law) is used to convert from power to voltage. The The duty cycle is the fraction of time that the transistor delivers power to the sensor. the Themainmaintask of task the of deck the driver deck reads driver the MOX reads sensor the MOX output voltage sensor and output temperature/ voltage and the humidity values The mainfrom temperature/humidity thevalues task SHT25 of and thefrom sends deckthe them SHT25 driver to the and reads theground sends themstation MOX at to the ground sensor 10 Hz. output station voltage atand 10 Hz. the temperature/humidity values from the SHT25 and sends them to the ground station at 10 Hz. Figure Figure 4.4.Schematic Schematic of the of the conditioningelectronic electroniccircuit circuit for each MOX sensorthe in MOX the MOX deck,deck, Figure 4. Schematic of theconditioning conditioning electronic circuit forfor each each MOXMOX sensor sensor in theinMOX deck, using using using PWM for powering and a voltage divider PWM for powering and a voltage divider for read-out. for read-out. PWM for powering and a voltage divider for read-out. 2.2. Experimental Arena, Gas Source and External Localization System 2.2. 2.2. Experimental Experimental Arena,Gas Arena, GasSource Source and External ExternalLocalization System Localization System All experiments were performed in a large robotics laboratory (160 m2 × 2.7 m height) at Örebro All All University experiments (Sweden). experiments were performed Theperformed were inisaa divided laboratoryin large largerobotics into laboratory robotics three (160(160 connected laboratory m2 ×areas m2.7 2 ×m height) (R1–R3) 2.7 atof m height)Örebro 132 m2 and at Örebro University a contiguous (Sweden). The laboratory 2 is divided into three connected areas (R1–R3) of 132 m 2 and a University room (R4) of (Sweden). 28 laboratory The m (Figureis5).divided The ventilation into threesystem of theareas connected laboratory (R1–R3)wasof not 132modified m2 and a contiguous room (R4) of 28 m22 (Figure 5). The ventilation system of the laboratory was not modified for the experiments and all windows and doors were kept closed. contiguous room (R4) of 28 m (Figure 5). The ventilation system of the laboratory was not modified for the experiments and all windows and doors were kept closed. for the experiments and all windows and doors were kept closed. Figure 5. Experimental arena. (left) Frontal picture; (right) Schematic top view. The green squares Figure 5. Experimental arena. (left) Frontal picture; (right) Schematic top view. The green squares indicate the position of the UWB anchors, which are positioned along two inverted triangles indicate the position of the UWB anchors, which are positioned along two inverted triangles (green (green lines). lines). Figure 5. Experimental arena. (left) Frontal picture; (right) Schematic top view. The green squares Toindicate obtain To obtaintheposition the 3D position the 3D of theofUWB position the of thedrone, drone,we anchors, we used which anexternal usedare an externallocalization positioned localization along two system inverted system (Loco triangles (Loco positioning (green positioning system, Bitcraze lines). system, AB) [66] Bitcraze AB) based on ultra-wide [66] based band (UWB) on ultra-wide radio transmitters. band (UWB) The system radio transmitters. is composed The system is of sixcomposed anchors of that six are positioned anchors in knowninlocations that are positioned of the room known locations and one of the room and tag thatthat one tag is fixed is fixedto the drone. The toTo anchors theobtain drone. the The were 3D placed position anchors inplaced of were thedrone, the central we in the area of the used central laboratory, an external area shaped localization of the laboratory, insystem shaped twoininverted (Loco two triangles positioning inverted triangles system, (below and (below Bitcraze above AB) and above the [66] based the flight flight area), as area), as on recommended recommended ultra-wide band by (UWB) by the the manufacturer. radio transmitters. manufacturer. The The tag ontag on Thethe thedroneis system composed of six anchors that are positioned in known locations of the room and one tag that is fixed to the drone. The anchors were placed in the central area of the laboratory, shaped in two inverted triangles (below and above the flight area), as recommended by the manufacturer. The tag on the Sensors 2019, 19, 478 8 of 25 Sensors 2018, 18, x FOR PEER REVIEW 8 of 26 continuously sends short high frequency radio messages to the anchors and estimates its relative position to them based drone continuously on the sends timestamps short of transmitted high frequency and received radio messages messages. to the anchorsThe andaccuracy estimates in the its estimated positiontois them relative position approximately based on10the cmtimestamps if the tag is within the space and of transmitted delimited receivedby the anchors and messages. The there accuracyis line-of-sight (LOS)position in the estimated betweenisthe anchors and10 approximately the cm tag [67].tag if the In is the absence within the of these space conditions, delimited by the system may still work but with degraded accuracy. the anchors and there is line-of-sight (LOS) between the anchors and the tag [67]. In the absence of theseAconditions, gas leak the wassystem emulatedmay by stillplacing work but a with smalldegraded beaker filled accuracy.with 200 mL of ethanol 96% (Sigma–Aldrich, A gas leak was emulated by placing a small beaker filled with 200 Ethanol Germany) in different locations of the arena (Figure 6). was used mL of ethanol 96%because (Sigma– it is non-toxic Aldrich, and easily Germany) detectable in different by MOX locations ofsensors. the arena Three (Figureexperiments 6). Ethanol were wascarried out to check used because the it is non- viability of the proposed system for GSL in complex environments. In the toxic and easily detectable by MOX sensors. Three experiments were carried out to check the viability first experiment, the gas source was placed of the proposed on top system forofGSL a table (height environments. in complex = 1 m) in the small In theroom first (R4). In the second experiment, the gasexperiment, source was the source placed wasofplaced on top a table inside the=suspended (height 1 m) in theceiling small (height room (R4). = 2.7Inm)thenear the entrance second experiment, to the thelab (R1). source Since the piping system of the lab runs through the suspended ceiling, the was placed inside the suspended ceiling (height = 2.7 m) near the entrance to the lab (R1). Since the gas source could represent a leak in piping one of system of the the pipes. lab runsInthrough these two the experiments, suspended ceiling, a 12 V theDCgasfan (Model: source couldAD0612HB-A70GL, represent a leak in ADDA one of the Corp., Taiwan) pipes. In theseplaced behind the beaker two experiments, a 12 V DC facilitated fan (Model:the dispersion AD0612HB-A70GL,of the chemicals in the ADDA Corp., environment, creating a plume. In the third experiment, the source was Taiwan) placed behind the beaker facilitated the dispersion of the chemicals in the environment, placed inside a power outlet box (height creating = 0.9 m) a plume. andthird In the a fishexperiment, tank bubbler thewas used source wasto increase the evaporation placed inside a power outlet rate.box The source= (height in 0.9this location m) and a fishcould simulatewas tank bubbler theused earlytostages increase of the an evaporation electrical firerate.(mostTheofsource them in arethis caused locationby faulty electrical outlets) where volatile organic compounds (VOCs) are released could simulate the early stages of an electrical fire (most of them are caused by faulty electrical into the environment. The threewhere outlets) experiments volatilestarted organicfive minutes after compounds setting (VOCs) areup the source released intoandthe turning on the DC environment. The fan or three the bubbler. experiments started five minutes after setting up the source and turning on the DC fan or the bubbler. Figure Figure 6. Gassource 6. Gas sourcelocation locationininthe thethree threeexperiments. experiments.(a) (a) Experiment Experiment 1: inside 1: inside small small room; room; (b) (b) Experiment 2: hidden in suspended ceiling; (c) Experiment 3: hidden in a power outlet box. Experiment 2: hidden in suspended ceiling; (c) Experiment 3: hidden in a power outlet box. 2.3. Gas Sensor Calibration and Limit of Detection (LOD) Estimation Sensors 2019, 19, 478 9 of 25 2.3. Gas Sensor Calibration and Limit of Detection (LOD) Estimation The MOX sensor was calibrated under laboratory conditions, to compensate the non-linear response and obtain measurements in concentration units. Ethanol concentrations up to 50 ppm were generated Sensors using 2018, 18,the permeation x FOR PEER REVIEWmethod [68], humidified to 30 % r.h. and delivered in random 9 of 26order at 70 mL/min to the chamber containing the sensor under test. The uncertainty of the concentration The MOX sensor was calibrated under laboratory conditions, to compensate the non-linear reaching the gas chamber was determined by propagation of the main sources of error, namely response and obtain measurements in concentration units. Ethanol concentrations up to 50 ppm were the permeation rate, the oven temperature and the MFCs. Uncertainty values ranged from 50 ppb generated using the permeation method [68], humidified to 30 % r.h. and delivered in random order (at 1 ppm at 70 of ethanol) mL/min to 1chamber to the ppm (atcontaining 50 ppm of ethanol), the which sensor under represents test. a relative The uncertainty uncertainty of the of 2–5%. concentration The concentration range and humidity level for calibration samples were selected reaching the gas chamber was determined by propagation of the main sources of error, namely thebased on previous experience by therate, permeation authors performing the oven similar temperature experiments and the in the same MFCs. Uncertainty test values environment. ranged from 50 ppb (at 1 To ppmestimate the LOD of ethanol) (ppm), to 1 ppm weppm (at 50 usedofthe simplified ethanol), whichLOD formula represents [69]: uncertainty of 2–5%. a relative The concentration range and humidity level for calibration samples were selected based on previous experience by the authors performing similar 3.3 × s0 in the same test environment. LODexperiments = , (2)  LOD formula [69]: To estimate the LOD (ppm), we used the simplified where s0 is the estimated standard deviation of LODblank3.3 𝑠 measurements , (assuming homoscedasticity (2) and 𝐴 normality) and  is the estimated slope of the calibration graph (assuming linearity). To account for where 𝑠 variability the expected is the estimated betweenstandard deviation of the calibration blank setup measurements and (assuming [70], the test environment homoscedasticity we estimated s0 from and normality) and a preliminary 𝐴 is the estimated exploration slope scenario of the target of the calibration graph (assuming in the absence linearity). of the gas. The LODTo account was used for the expected variability between the calibration setup and the test environment [70], during the gas source localization experiments to remove false alarms, by setting to zero any measuredwe estimated 𝑠 from abelow concentration preliminary exploration of the target scenario in the absence of the gas. The LOD was used the LOD. during the gas source localization experiments to remove false alarms, by setting to zero any measuredofconcentration 2.4. Detection ‘Bouts’ below the LOD. To 2.4.compute Detection of the ‘bouts’ from the MOX response, we adapted the signal processing pipeline ‘Bouts’ proposed To bycompute Schmuker the et al. [65]. ‘bouts’ fromThe the goal MOXof this algorithm response, we adapted is totheextract signal the rising edges processing pipeline of the smoothed derivative proposed of the sensor by Schmuker response, et al. [65]. The goalwhich are called of this algorithm‘bouts’. is to Schmuker’s algorithm extract the rising edges isofbased the on a non-causal smoothed Gaussian derivative smoothing of the sensorfilter that prevents response, which arereal-time bout detection, called ‘bouts’. Schmuker’s includes algorithman unnecessary is based derivative, embeds the on a non-causal computation Gaussian smoothing of the derivative filter within that prevents the smoothing real-time filter (which bout detection, includesmay an lead unnecessary derivative, embeds the computation of the derivative within to potential implementation errors) and does not smooth the second derivative of the signal where the smoothing filter (which may lead the ‘bouts’ areto potential implementation segmented. The derivative errors) that and does not converts s tosmooth the is x in [65] second derivative unnecessary of the signal because the emaα where the ‘bouts’ are segmented. The derivative that converts 𝑠 to transformation already differentiates the input signal. The source code published by Schmuker 𝑥 in [65] is unnecessary because et al. is the 𝑒𝑚𝑎 transformation already differentiates the input signal. The source code published by not affected by this error because the transformation is implemented by calling the Python function Schmuker et al. is not affected by this error because the transformation is implemented by calling the pandas.ewma(), which provides the functionality of an EWMA filter (i.e. it does not perform the Python function pandas.ewma(), which provides the functionality of an EWMA filter (i.e. it does not derivative). perform We the replaced derivative).the WeGaussian replaced smoothing the Gaussianfilter with a filter smoothing causal (realizable) with exponentially a causal (realizable) weighted exponentially weighted moving average (EWMA) filter, removed the unnecessary derivative, the moving average (EWMA) filter, removed the unnecessary derivative, decoupled derivative fromthe decoupled thederivative EWMA from filter the andEWMAsmoothed the second filter and smoothed derivative. the second We chose the derivative. We EWMA chose thefilter because EWMA it is filter causal, easyitto because implement is causal, easy toand is the same implement and isfilter used the same by used filter Schmuker et al. ettoal.smooth by Schmuker to smooth thein the derivative derivative their boutin their bout computation computation algorithm. algorithm. TheThe proposed proposed boutcomputation bout computation pipeline pipeline is is represented represented in Figure in Figure 7. 7. Figure 7. Flow Figure diagram 7. Flow of the diagram improved of the improvedbout boutcomputation. Themeaning computation. The meaningof of each each symbol symbol is given is given in in the text. the text. The sensor response The sensor x is first response 𝑥 issmoothed usingusing first smoothed a EWMA low-pass a EWMA filter to low-pass remove filter high-frequency to remove high- At time t,noise. noise.frequency At time 𝑡, the the smoothed value xs (t) isvalue smoothed found 𝑡 computing 𝑥 by is found by computing 𝑥 𝑡 1 𝛼 ∙𝑥 𝑡 1 +𝛼∙𝑥 𝑡 , (3) xs (t) = (1 − α)· xs (t − 1) + α· x (t), (3) where 𝑥 𝑡 is the observation at time t, 𝑥 𝑡 1 is the previous output of the filter and the smoothing factor 𝛼 (0 < 𝛼 ≤ 1) controls the speed at which older responses are dampened. The Sensors 2019, 19, 478 10 of 25 where x (t) is the observation at time t, xs (t − 1) is the previous output of the filter and the smoothing factor α (0 < α ≤ 1) controls the speed at which older responses are dampened. The smoothing factor α in the EWMA filter is equivalent to the cut-off frequency of a low-pass filter. For example, Pashami et al. [71] found out that α = 0.9 represents a cut-off frequency of 0.44 Hz. It is convenient to express the smoothing factor as a function of the half-life time τhalf (s) ! log 0.5 α = 1 − exp , (4) τhal f · f s where f s is the sampling frequency (Hz) of x. The smoothed response is differentiated and smoothed again to increase the signal-to-noise ratio (SNR), producing xs0 . The ‘bouts’ are the rising edges of xs0 , 00 which are delimited by two consecutive zero-crossings in the positive derivative10of Sensors 2018, 18, x FOR PEER REVIEW x 0 , i.e., xs > 0. of 26 s The amplitude of a bout is defined as xs0 at the end of the respective bout segment minus xs0 at the smoothing factor 𝛼 in the EWMA filter is equivalent to the cut-off frequency of a low-pass filter. For start of the same bout segment. To remove low-amplitude bouts produced by noise, Schmuker et al. example, Pashami et al. [71] found out that 𝛼 = 0.9 represents a cut-off frequency of 0.44 Hz. It is propose to estimate convenientthetonoise expressthreshold (bfactor the smoothing thr ) using the 3-sigma as a function rule of the half-life time 𝜏 (s) . 𝛼 1 exp , (4) bthr = µ + 3σ,∙ (5) where 𝑓 is the sampling frequency (Hz) of 𝑥 . The smoothed response is differentiated and where µ and σsmoothed are theagain to increasemean estimated the signal-to-noise and standardratio (SNR), 𝑥 . The ‘bouts’ areof producing respectively, deviation, thethe rising distribution of edges of 𝑥 , which are delimited by two consecutive zero-crossings in the positive derivative of 𝑥 , amplitudes of i.e., bouts𝑥 detected in the sensor baseline (i.e., in the absence of gas). Bouts 0. The amplitude of a bout is defined as 𝑥 at the end of the respective bout segment with amplitude minus 𝑥 lower than bthr are filtered out. The algorithm that we propose reduces the number of parameters from at the start of the same bout segment. To remove low-amplitude bouts produced by noise, Schmuker et al. propose to estimate the noise threshold (𝑏 ) using the 3-sigma rule three to two (τhalf and bthr ), which we estimate after a preliminary exploration of the target scenario in 𝑏 𝜇 + 3𝜎, (5) the absence of gas (i.e., using the signals corresponding to blank measurements). where μ and σ are the estimated mean and standard deviation, respectively, of the distribution of amplitudes of bouts detected in the sensor baseline (i.e., in the absence of gas). Bouts with amplitude 2.5. Effect of Rotors on MOX Sensor Signals lower than 𝑏 are filtered out. The algorithm that we propose reduces the number of parameters from three to two (𝜏 and 𝑏 ), which we estimate after a preliminary exploration of the target Previous work using the CF2 indicate that turbulence generated by the propellers may severely scenario in the absence of gas (i.e., using the signals corresponding to blank measurements). affect the MOX sensor signals [29]. To evaluate this effect in our platform, we performed a set of measurements2.5. Effectaof gas near Rotorssource on MOX under Sensor Signals two conditions: rotors switched on or rotors switched off. The drone was placed Previous onwork using the CF2 indicate that a height-adjustable turbulence generatedto stand—designed byminimize the propellersits may severely interference with the affect the MOX sensor signals [29]. To evaluate this effect in our platform, we performed a set of rotors’ airflow—that could be moved around the source (Figure 8a). Using a stand is necessary to measurements near a gas source under two conditions: rotors switched on or rotors switched off. The perform measurements drone was placedwhen on the rotors are switched a height-adjustable off, to stand—designed and we also minimize used it to its interference perform with the rotors’measurements with the rotorsairflow—that switched could on. Thebe moved around the source (Figure 8a). Using a stand is necessary to perform gas source was an open ethanol bottle (Figure 8b) and measurements measurements when the rotors are switched off, and we also used it to perform measurements with were performed theabove the source rotors switched on. The (vertical gas sourcedistance was an openbetween 25 and ethanol bottle 65 8b) (Figure cm)andand in front of it (at 50 cm) measurements At each location, werethe sensorabove performed response was(vertical the source recorded forbetween distance 25 min 25 (using and 65 cm) anand external in front battery), of it (at 50 first with the cm) At each location, the sensor response was recorded for 25 min (using an external battery), first rotors switched off and, after cleaning the room, measurements were repeated with the rotors switched off and, after cleaning the room, measurements were repeated with the four with the four rotors spinning at 10,000 rotorsrpm spinning(this is typical at 10,000 for rpm (this hovering). is typical for hovering). Figure 8. Setup for 8. Figure assessing the effect Setup for assessing of the the effect rotors of the onthethe rotors on MOX MOX sensorsensor signals. signals. (a) Top view(a) Top view of the of the stand used to hold the drone at different heights while minimizing interference with the rotors air stand used to hold the drone at different heights while minimizing interference with the rotors air flow; flow; (b) Photo of an experiment with the drone placed 25 cm above an ethanol bottle (gas source), (b) Photo of anoverlaid experiment with theofdrone with an illustration placed 25 cm above an ethanol bottle (gas source), overlaid a gas cloud. with an illustration of a gas cloud. Sensors 2019, Sensors 19,18, 2018, 478x FOR PEER REVIEW 1111ofof2526 2.6. Gas Source Localization Strategies 2.6. Gas Source Localization Strategies Two gas source localization strategies, one based on the instantaneous response and the other Two gas source localization strategies, one based on the instantaneous response and the other one based on odor hits are evaluated using the nano-drone. In both cases, the drone was sent to fly one based on odor hits are evaluated using the nano-drone. In both cases, the drone was sent to fly along a predefined sweeping path consisting of two 2D rectangular sweepings at different heights along a predefined sweeping path consisting of two 2D rectangular sweepings at different heights (0.9 m and 1.8 m), collecting measurements in motion (Figure 9). These two heights divide the vertical (0.9 m and 1.8 m), collecting measurements in motion (Figure 9). These two heights divide the vertical space of the lab in three parts of equal size. Flying first at a lower altitude minimizes the impact of space of the lab in three parts of equal size. Flying first at a lower altitude minimizes the impact of the the propellers’ downwash in the gas distribution. For safety reasons, the trajectory is designed to propellers’ downwash in the gas distribution. For safety reasons, the trajectory is designed to ensure ensure enough clearance around obstacles and walls, and people working inside the laboratory were enough clearance around obstacles and walls, and people working inside the laboratory were told to told to remain in their seats during the experiments. The ground station communicates the flight path remain in their seats during the experiments. The ground station communicates the flight path to the to the drone as a sequence of (𝑥, 𝑦, 𝑧) waypoints, with a target flight speed of 1.0 m/s. drone as a sequence of (x, y, z) waypoints, with a target flight speed of 1.0 m/s. As the drone navigates the environment, it reports the instantaneous concentration and its As the drone navigates the environment, it reports the instantaneous concentration and its location to the ground station. At the end of the exploration, the ground station uses all the received location to the ground station. At the end of the exploration, the ground station uses all the received information to compute a 3D map of the instantaneous response (first strategy) and the bout information to compute a 3D map of the instantaneous response (first strategy) and the bout frequency frequency (second strategy). The location of the gas source is estimated in both cases as the cell of the (second strategy). The location of the gas source is estimated in both cases as the cell of the map with map with maximum value. We will also discuss the viability of both methods for real-time plume maximum value. We will also discuss the viability of both methods for real-time plume tracking, tracking, assuming that the drone would follow the gradient of instantaneous concentration or the assuming that the drone would follow the gradient of instantaneous concentration or the gradient of gradient of bout frequency (computed using a sliding window of 5 s). bout frequency (computed using a sliding window of 5 s). Figure 9. 9.Predefined Figure Predefinednavigation navigationstrategy strategybased basedononzig-zag zig-zagsweeping sweepingatattwo twoheights heights(0.9 (0.9and and1.8 1.8m). m). The green squares indicate the location of the UWB anchors. The green squares indicate the location of the UWB anchors. 3. Results 3. Results 3.1. Calibration, LOD and Optimum Parameters for Bout Detection 3.1. Calibration, LOD and Optimum Parameters for Bout Detection A preliminary exploration of the test environment in the absence of gas was performed to estimate A preliminary the sensor noise and theexploration optimumof theparameters bout test environment (Figurein the The 10A). absence of gas was raw response wasperformed smoothedto estimate the sensor noise and the optimum bout parameters (Figure 10A). using an EWMA filter with a smoothing factor τhalf = 0.25 s. The noise can be approximated byThe raw response was smoothed using an EWMA filter with a smoothing − 1 factor 𝜏 0.25 a Gaussian distribution with mean value of 2.047 MΩ and standard deviation of 0.013 MΩ 1be s . The noise can− approximated (Figure 10B). Theby a Gaussian observed distribution variability with in was used mean value of 2.047 combination MΩcalibration with the and standard deviation line (Figure 11A)of 0.013 MΩ (Figure 10B). The observed variability was used in combination with the to estimate the LOD (Equation (2)), assuming homoscedasticity. The calibration line behaved linearly calibration line (Figure 11A) to estimate the LOD (Equation 2), assuming homoscedasticity. in the range 1–50 ppm after applying the logarithm to both concentration and response. The bout The calibration line behaved threshold amplitude linearly in(bthe) was range 1–50 ppm estimated afterppm/s as 0.04 applying the logarithm by applying to both Equation (5) to concentration the amplitude and of thr response. The bout amplitude threshold (𝑏 ) was the bouts detected in the calibrated blank signals (Figure 11B).estimated as 0.04 ppm/s by applying Equation 5 to the amplitude of the bouts detected in the calibrated blank signals (Figure 11B). Sensors 2018,19, Sensors2019, 18,478 x FOR PEER REVIEW 12 of 12 of 25 26 Sensors 2018, 18, x FOR PEER REVIEW 12 of 26 Figure 10. (A) 2D map of MOX sensor response during 15 minutes of random exploration of the target map of MOX Figure 10. (A) 2D map MOX sensor response during 15 min of random minutes exploration of random of theoftarget exploration area the target area without gas; (B) Histogram of blank readings, with a Gaussian curve Ν(𝜇, 𝜎 ) superimposed. 2 ) superimposed. without gas; (B) area without gas;Histogram of blank (B) Histogram readings, of blank with awith readings, Gaussian N(µ, σΝ(𝜇, curvecurve a Gaussian 𝜎 ) superimposed. Figure Figure11.11.(A) Calibration (A) line in Calibration thein line range the 1–50 ppm range (log-log 1–50 ppmplot), with blank (log-log plot),variability with blank superimposed variability Figure 11. (A) Calibration line in the range 1–50 ppm (log-log plot), with blank variability at each concentration superimposed at eachlevel (see inset).level concentration The (see LODinset). is estimated The LOD using Equation using is estimated (2); (B)Equation Histogram of 2; (B) superimposed at each concentration level (see inset). The LOD is estimated using Equation 2; (B) amplitudes Histogram of bouts detected in the calibrated blank signals. amplitudes of bouts detected in the calibrated b thr blank signals. 𝑏 (Equation (5)) is indicated by a red (Equation 5) is Histogram of amplitudes of bouts detected in the calibrated blank signals. 𝑏 (Equation 5) is dashed vertical indicated line.dashed vertical line. by a red indicated by a red dashed vertical line. 3.2. Effect of Propulsion on MOX Signals 3.2. Effect of Propulsion on MOX Signals 3.2. Effect of Propulsion on MOX Signals When the propellers were switched off (Figure 12A), the concentration fluctuations due to the gas When the propellers were switched off (Figure 12A), the concentration fluctuations due to the When the evaporating frompropellers the ethanol were switched bottle off (Figure were clearly 12A), reflected concentration on the on-board sensorfluctuations signals at 25due to the cm above gas evaporating from the ethanol bottle were clearly reflected on the on-board sensor signals at 25 gas source the evaporating (greenfromtrace)the ethanol and, to lessbottle extent,were at 50clearly cm inreflected front of the on the source on-board sensor and (blue signals trace). At 65atcm25 cm above the source (green trace) and, to less extent, at 50 cm in front of the source and (blue trace). cm above above the source (green trace)theand, to less extent, not at 50 cm in fronttheof source the source andimmediately (blue trace). At 65 thecm source above (yellow the source trace), (yellowsensor trace), does the sensor seem doesto detect not seem to detect except the source except At 65opening after cm above the source the bottle (t = 2(yellow trace), min). When thethe samesensor does not seem measurements to detect with were repeated the source except the propellers immediately after opening the bottle (𝑡 2 min). When the same measurements were repeated with immediately switched after opening on (Figure 12B), onthe bottle (𝑡 2 min). thefluctuations Whenatthe 25ofsame measurements were repeated with the propellers switched (Figure 12B),ofthe thefluctuations signals cmtheabove signalsthe at source 25 cm became abovelesstheintense source the propellers but moreless switched on frequent. (Figure 12B),observed the fluctuations of13,thewhere signalsthe at bouts 25 cmdetected above the source became intense This but morecan be better frequent. This can in be Figure better observed in Figure 13, where the with the bouts became less propellers intense on switched butaremoreas frequent. twice as This can frequent be better than when observed the in Figure propellers are 13, where switched off the bouts (Table 1). detected with the propellers switched on are as twice as frequent than when the propellers are detected This with effectoff the is (Table propellers even more switched noticeable on are as twice as frequent than when the propellers are switched 1). This effect is at 65 cm even moreabove the source noticeable at 65because cm above the bout the frequency source becauseincreases the bout switched from 0.48 off to (Table 7.74 1). This effect bouts/min by is even more switching on noticeable the at 65improving propellers, cm above the the source detection because of thethe bout source. frequency increases from 0.48 to 7.74 bouts/min by switching on the propellers, improving the frequency The increases propellers produced fromnonetheless 0.48 to 7.74a negative bouts/min by switching effect in the on acquired signals the propellers, in front improving of the the source, detection of the source. The propellers produced nonetheless a negative effect in the signals acquired detection the reducing of the source. The propellers produced two. nonetheless a negative effect in the signals acquired in front of thebout frequency source, reducing by athe factor boutoffrequency by a factor of two. in front of the source, reducing the bout frequency by a factor of two. Sensors 2019, 19, 478 13 of 25 Sensors 2018, 18, x FOR PEER REVIEW 13 of 26 Sensors 2018, 18, x FOR PEER REVIEW 13 of 26 Figure 12. Figure Sensor Sensorsignals signals(log (logscale) near scale) an an near evaporating evaporating source. (A) (A) source. Propellers switched Propellers off; off; switched (B) Figure 12. 12. Sensor signals (log scale) near an evaporating source. (A) Propellers switched off; (B) Propellers (B) switched Propellers on. switched The ethanol bottle is opened at t = 2 min. Propellers switched on. on. TheThe ethanol ethanol bottle bottle is opened is opened at tat = 2t = 2 min. min. 0 ′ Figure 7) of the sensor signals at 50 cm in front of the source Figure 13. Smoothed Figure 13. Smoothedderivative (i.e.,xs𝑥in derivative(i.e., in Figure 7) of the sensor signals at 50 cm in front of the Figure 13. Smoothed derivative (i.e., 𝑥𝑠′𝑠 in Figure 7) of the sensor signals at 50 cm in front of the (blue line), 65 cm above the source (yellow source (blue line), 65 cm above the source (yellow line) and 25and line) cm 25 above the source cm above (green(green the source line). line). BoutsBouts with source (blue line), 65 cm above the source (yellow line) and 25 cm above the source (green line). Bouts amplitude higher than + 3σ are highlighted in red. In the left column, the propellers with amplitude higher than 𝜇 + 3𝜎 are highlighted in red. In the left column, the propellers are µ are switched off with amplitude higher than 𝜇 + 3𝜎 are highlighted in red. In the left column, the propellers are whereas in the right column they are switched on. The ethanol bottle is opened at t = switched off whereas in the right column they are switched on. The ethanol bottle is opened at 𝑡 2 2 min. switched off whereas in the right column they are switched on. The ethanol bottle is opened at 𝑡 2 min. 1. Characterization of MOX signals at different distances of the source under two conditions: Table min. propellers switched on or off. Table 1. Characterization of MOX signals at different distances of the source under two Table 1. Characterization of MOX signals at different distances of the source under two conditions: Variance Bout Frequency Bout Amplitude Distance propellers conditions: propellers switched switched Propellers on or on or Mean off. off. (ppm) (ppm2 ) (Bouts/min) (ppm/s) Bout Frequency Bout Amplitude Distance Propellers OFF Mean (ppm) 10.05 Variance (ppm2) Bout3.52 60.46 Frequency Bout Amplitude 0.39 Distance Above 25 cm Propellers Mean (ppm) Variance (ppm2) (Bouts/min) (ppm/s) ON 9.22 29.97 7.69 (Bouts/min) 0.084 (ppm/s) OFF OFF 10.051.39 60.46 0.053 3.52 0.48 0.39 0.027 Above OFF 10.05 60.46 3.52 0.39 Above25 Above 65cm 25 cm cm ON 9.22 29.97 7.69 0.084 ON ON 9.222.67 0.53 29.97 7.74 7.69 0.015 0.084 OFF OFF 1.391.68 0.053 0.59 0.48 1.13 0.027 0.10 Above 65 cm Front 50 cm OFF 1.39 0.053 0.48 0.027 Above 65 cm ON ON 2.671.45 0.53 0.12 7.74 0.47 0.015 0.10 ON 2.67 0.53 7.74 0.015 OFF 1.68 0.59 1.13 0.10 Front 50 cm OFF 1.68 0.59 1.13 0.10 Front 50 cm ON 1.45 0.12 0.47 0.10 ONof the average The steady increase 1.45 concentration 0.12and higher bout 0.47 0.10 when the frequency observed drone was sampling above the source with the rotors switched on can be explained based on how The steady increase of the average concentration and higher bout frequency observed when the The steadyinteract the propellers increase of the with theaverage concentration gas surrounding the and higher drone. Thebout frequency propellers observed generate when the a downwash drone was sampling above the source with the rotors switched on can be explained based on how drone was sampling above the source with the rotors switched on can be explained based on how (i.e., a downward airflow) that acts as an opposing force to the gas moving upwards by convection, the propellers interact with the gas surrounding the drone. The propellers generate a downwash (i.e., the propellers interact with the gas surrounding the drone. The propellers generate a downwash (i.e., a downward airflow) that acts as an opposing force to the gas moving upwards by convection, a downward airflow) that acts as an opposing force to the gas moving upwards by convection, Sensors 2019, 19, 478 14 of 25 Sensors 2018, 18, x FOR PEER REVIEW 14 of 26 breaking the breaking the laminar flow flowinto intoaaturbulent turbulentgasgascloud that cloud spreads that around spreads around the the drone’s frame drone’s (Figure frame 14A), (Figure probably increasing the background concentration. When the gas 14A), probably increasing the background concentration. When the gas is spread is spread around the drone, the drone, thepropellers the propellersdrag dragthe thegas gaspatches patchestowards towardsthe theMOX MOXsensor sensor(Figure (Figure14B)14B)and andthis thisprobably probablyincreases increases thebout the boutfrequency. frequency. Therefore, Therefore,from fromthis thisstudy studywe wecan canconclude concludethat thatthethepropellers propellersseverely severelychange change the gas distribution near the drone, which in turn affects the MOX signals. However, the gas distribution near the drone, which in turn affects the MOX signals. However, it may still it may stillbe be possibleto possible toextract extractrelevant relevant features features for for gas gas source source localization localization (e.g., (e.g., the the bout bout frequency). Figure Figure 14. 14. Aerodynamics AerodynamicsofofCrazyflie Crazyflie2.02.0when when thethe four rotors four areare rotors spinning, spinning, visualized visualizedusing a using Deskbreeze wind tunnel (Courtesy of Bitcraze AB). The drone is fixed to one of the walls a Deskbreeze wind tunnel (Courtesy of Bitcraze AB). The drone is fixed to one of the walls of the of the tunnel using a 3D tunnel printed using stand and a 3D printed dryand stand ice fog dryisice emitted from (A)from fog is emitted below(A)the drone below theordrone (B) above or (B)the drone. above the Itdrone. showsItthe downwash shows of the propellers the downwash and how of the propellers and part howofpart the fog reaches of the the MOX fog reaches gas sensor the MOX (red gas sensor arrow). The MOX (red arrow). deck has The MOX deckbeen overlaid has been to thetooriginal overlaid images the original for visual images clarity. for visual clarity. 3.3.Experiment 3.3. Experiment1:1:Localization LocalizationofofaaSource Source17 17mmAway Awayfrom fromthe theStarting StartingPoint Point Inthe In thefirst firstexperiment, experiment,the thedrone dronetook tookoffoffnear nearthe theentrance entranceof ofthe thelab lab(R1), (R1),1717meters metersdownwind downwind ofaagas of gassource sourcelocated locatedin in the the other other end end of of the the laboratory laboratory (R4). (R4). Flying Flyingalong alongaapredefined predefinedtrajectory trajectory that took 2.75 min, the drone acquired measurements to build a map that took 2.75 minutes, the drone acquired measurements to build a map of the instantaneousof the instantaneous response (Figure 15A). response (FigureFrom thisFrom 15A). map this it is map evident it isthat the gas evident source that must the gas be inmust source R4, because be in R4,the maximum because the concentration (35 ppm) was found there while concentrations below maximum concentration (35 ppm) was found there while concentrations below 5 ppm were 5 ppm were measured in the rest of the lab. The gas plume can be outlined from the obtained bout map, which measured in the rest of the lab. The gas plume can be outlined from the obtained bout map, which shows the highest bout density shows the(25highest bouts/min) also in R4. bout density Most bouts also (25 bouts/min) wereindetected R4. Mostduring bouts the werefirst part ofduring detected the exploration the first (Figure 15C), when the drone was flying at the same height of the source part of the exploration (Figure 15C), when the drone was flying at the same height of the source(Figure 15B). Nonetheless, multiple (Figure bouts 15B). were also detected Nonetheless, multiple flying above bouts werethe source also detected (t = 110–120 in R4flying above thes). The cellsincorresponding source R4 (𝑡 110– to the 120 s). maximum concentrationto The cells corresponding andthemaximum maximumbout frequency and concentration weremaximum found at 0.94boutand 1.16 m of frequency the were true source location, respectively. found at 0.94 and 1.16 m of the true source location, respectively. Sensors 2019, 19, 478 15 of 25 Sensors 2018, 18, x FOR PEER REVIEW 15 of 26 Figure Figure 15. Resultsof 15. Results ofExperiment Experiment1.1.(A)(A) 2D2D mapmap of the of the instantaneous instantaneous concentration concentration (ppm) (ppm) in loginscale, log scale, with bouts represented by with bouts represented by blue circles ( 𝑏 blue circles (b thr = 0.52 ppm/s). A hand-drawn ellipse 0.52 ppm/s). A hand-drawn ellipse outlines the outlines the approximate approximate plumeplumeshapeshapebasedbased on location on the the location of bouts. of bouts. The average The average bout frequency bout frequency along along the y- the y-axis is shown in the leftmost panel. The box plot below the map axis is shown in the leftmost panel. The box plot below the map represents the instantaneousrepresents the instantaneous concentration concentrationalong alongthe thex-axis; x-axis;(B) (B)Trajectory Trajectoryof ofthe thedrone dronealong alongthe thez-axis. z-axis. (C) (C) Temporal Temporalevolution evolutionof of the instantaneous concentration (ppm) on a log scale, with detected bouts highlighted the instantaneous concentration (ppm) on a log scale, with detected bouts highlighted in red (the black in red (the black star starindicates indicatesthe thestart startofofa about). bout).The Theidentifiers identifiers R1–R4 R1–R4 between between panels panels (B)(B) and (C)(C) and indicate indicate thethe areaarea of the map in which the drone is flying at each moment. The maximum instantaneous of the map in which the drone is flying at each moment. The maximum instantaneous concentration concentration and the andmaximum the maximumbout bout frequency are indicated frequency by a green are indicated star and by a green staraand bluea triangle, respectively. blue triangle, respectively. For real-time plume tracking, following the gradient of instantaneous concentration would not For real-time plume tracking, following the gradient of instantaneous concentration would not easily lead the drone towards the source, since the concentrations measured in R1 are only slightly easily lead the drone towards the source, since the concentrations measured in R1 are only slightly higher than the LOD of the sensor (1.1 ppm). Although it is true that the instantaneous concentration higher than the LOD of the sensor (1.1 ppm). Although it is true that the instantaneous concentration increases when the drone crosses the plume (e.g., bouts #1, #2 and #6), it is not clear when the drone increases when the drone crosses the plume (e.g., bouts #1, #2 and #6), it is not clear when the drone exits the plume due to the increase in background concentration and slow sensor recovery. In this exits the plume due to the increase in background concentration and slow sensor recovery. In this sense, the bout frequency behaves like a more sensitive version of the instant response, exhibiting sense, the bout frequency behaves like a more sensitive version of the instant response, exhibiting abrupt changes when the drone enters/exits the plume. This can be seen at t = 85 s, when the drone abrupt changes when the drone enters/exits the plume. This can be seen at 𝑡 85 s, when the drone exits R4 and enters R3. The bout frequency sharply drops to zero, but the instantaneous concentration exits R4 and enters R3. The bout frequency sharply drops to zero, but the instantaneous concentration remains flat at 5 ppm. Another example is found by comparing both features when the drone crosses remains flat at 5 ppm. Another example is found by comparing both features when the drone crosses the door that connects R3 to R4 flying at h = 0.9 m (t = 75 s) or at h = 1.8 m (t = 110 s). In the first case, the door that connects R3 to R4 flying at ℎ 0.9 m (𝑡 75 s) or at ℎ 1.8 m (𝑡 110 s). In the the instantaneous concentration jumps from 4 to 35 ppm whereas in the latter case the step is only from first case, the instantaneous concentration jumps from 4 to 35 ppm whereas in the latter case the step is only from 5 to 10 ppm. On the other hand, the bout frequency shows a similar increase (0 to 23–25 Sensors 2019, 19, 478 16 of 25 5 toSensors 10 ppm.2018, On 18, xthe FORother PEER REVIEW hand, the bout frequency shows a similar increase (0 to 23–25 bouts/min) 16 of 26 in both cases. This could mean that the bout frequency, which simply counts events regardless of their bouts/min) in both cases. This could mean that the bout frequency, which simply counts events amplitude is more robust for detecting the source when its height is unknown. regardless of their amplitude is more robust for detecting the source when its height is unknown. To obtain the bout map shown in Figure 15A, we had to increase the bout threshold (bthr ) To obtain the bout map shown in Figure 15A, we had to increase the bout threshold (𝑏 ) from from 0.04 0.04 ppm/s ppm/s (computed (computed using using Equation Equation 5) to(5)) toppm/s 0.52 0.52 ppm/s (determined (determined from visual from visual inspection). inspection). The Theeffect effectofofthe the bout bout amplitude amplitude threshold threshold on on thethe bout bout mapmap is show is show in Figure in Figure 16. If16. theIfthreshold the threshold is set is set too toolowlow(blue (bluecircles), circles), slight variations of the background concentration can trigger slight variations of the background concentration can trigger bout detections that bout detections thatarearenotnot necessarily necessarily produced produced bygas by the theplume, gas plume, resulting resulting in adistribution in a sparse sparse distribution of bouts of bouts that hinders that hinders the localization the localization of theof source the source (localization (localization errorerror of 4.32 of 4.32 m).m). This This variability variability in inthethe background background concentration concentration maymay bebecaused causedby bythe thedownwash downwash of of the the drone, which whichbreaks breaksthe theplume plumeand anddisperses disperses fragments fragments of of thetheplume plume ininthe theroom. room.OnOnthethe other other hand, if if the the threshold thresholdisisset settoo toohigh high (green (green circles) circles) it isitvirtually is virtually guaranteedthat guaranteed thatevery everydetected detected bout bout is produced produced by bythe theplume. plume.However, However, ininthis case, this case, boutsbouts appear appear onlyonly nearthe near thesource, source,which whichreduces reduces the distance distanceat atwhich whichthe thesource sourcecan canbebe detected. detected. Only Onlywhenwhen thethe thresholdisissuitably threshold suitablyestimated estimated (Figure (Figure 15A), 15A), thethe bout bout map mapresembles resemblesananelongated elongated plumeplume thatthat cancanbebe used used totomove movetowards towardsthethesource source from from aa long long distance. distance.As Asofofnow, now,we weare are not aware not aware of a systematic approach to properly set this threshold without prior of a systematic approach to properly set this threshold without prior information about the sourceinformation about the source intensity. intensity. OneOne idea idea could could bebetotouse useaamulti-threshold multi-threshold approach, approach, i.e., i.e.,building buildingboutboutmapsmapswith with different different thresholds and searching for consistent source predictions over thresholds and searching for consistent source predictions over multiple threshold values.multiple threshold values. Figure Figure Effect 16.16. Effectofofthe thebout boutamplitude amplitudethreshold threshold inin the the results Experiment1.1.The results of Experiment Theblue bluecircles circles represent representbouts boutswithwithamplitude amplitudehigher higherthan than 0.04 0.04 ppm/s ppm/s (𝜇 (µ++3𝜎 3σ threshold) threshold)andandthe thegreen greencircles circles represent represent bouts bouts withwith amplitude amplitude higher higher than than 1.01.0 ppm/s. ppm/s. InIn each each case, case, a hand-drawn a hand-drawn ellipse ellipse outlinesthe outlines the approximate plume shape based on the location of the bouts. The green approximate plume shape based on the location of the bouts. The green and blue stars indicate and blue stars indicatethe the source source locationlocation estimate estimate in case, in each each case, according according to thetomaximum the maximum boutbout frequency. frequency. 3.4.3.4. Experiment 2: 2: Experiment Localization Localizationofofa aSource SourceHidden Hidden in in the the Suspended Ceiling(h(h==2.7 Suspended Ceiling 2.7m)m) In this experiment, In this experiment, the thesource sourcewas waslocated locatedjust just above the starting above the startingpoint pointofofthe theexploration, exploration, hidden hidden in the suspended ceiling (Figure 17). The resulting maximum concentration in the suspended ceiling (Figure 17). The resulting maximum concentration in the test room was in the test room was measured when the drone flew at ℎ 1.8 m, highlighting the importance of sampling in 3D forfor measured when the drone flew at h = 1.8 m, highlighting the importance of sampling in 3D localization localization andandmapping mappingofofelevated elevatedgas gassources. However, since sources. However, sincethe thesource sourceisispresumably presumablynot not directly directly exposed exposed toto thethe environment,concentrations environment, concentrations below below 33 ppmppm were werefound foundininmost mostlocations locationsof of thethe room, room, which which complicatesthe complicates theGSL GSLtask. task. The 𝜇+ The µ + 3𝜎 threshold yielded 3σ threshold yieldedbetter betterresults resultsthan thaninin thethe previous previous experiment, experiment, outlininga afeasible outlining feasibleplume plume shape shape according according to tothe thetrue truegas gassource sourcelocation locationand and predominant predominant windwind direction direction (Figure (Figure 18A).18A). The The instantaneous instantaneous response response map map andbouts and the the bouts map map suggest thatsuggest the gasthat the gas source is insource is in thebetween the division divisionR1 between and R2, R1which and R2, which represents represents a localization a localization error error of 4.0 and of 4.0 and 3.31 m, respectively. The localization error can be reduced to 2.22 m by 3.31 m, respectively. The localization error can be reduced to 2.22 m by increasing the bout threshold increasing the bout threshold from 0.04 ppm/sfromto0.04 0.18ppm/s ppm/s to (Figure 0.18 ppm/s (Figure 19). In 19). the this case, In this lowercase, the lower inter-bout inter-bout interval betweeninterval bouts between bouts #3 and #4 (as compared to bouts #1 and #2) leads to a higher bout frequency near the #3 and #4 (as compared to bouts #1 and #2) leads to a higher bout frequency near the source. source. Sensors2019, Sensors 2018,19, 18,478 x FOR PEER REVIEW 17ofof2526 17 Sensors 2018, 18, x FOR PEER REVIEW 17 of 26 Figure 17. 3D map of the instantaneous concentration (ppm) in Experiment 2. The black square Figure 17.17.3D3D mapofofthe map theinstantaneous instantaneous concentration concentration (ppm) in in Experiment 2.2.The black blacksquare Figure indicates the gas source location (x,y,z) = (14.0, 5.2, 2.7) m,(ppm) the black Experiment arrow the wind The direction square (positive indicates indicates the the gas gas source source location location(x,y,z) = (x,y,z) (14.0, = 5.2, (14.0, 2.7) 5.2, m, the 2.7) black m, thearrow blackthe wind arrow direction the wind (positive direction x-axis) and the letter x-axis) x-axis) and theand letter ‘S’ the starting point pointof ofthe thedrone drone(x,y,z) (x,y,z)= =(13.5, (13.5, 5.2, 0.0) m. (positive the‘S’letter the starting ‘S’ the starting point of the drone 5.2, (x,y,z) = 0.0) m.5.2, 0.0) m. (13.5, Figure 18. Results of Experiment 2. (A) 2D map of the instantaneous concentration (ppm), with odor Figurerepresented hits Figure 18.Results 18. by Resultsof blue circles (𝑏 ofExperiment Experiment 2.2.(A) 0.04 (A)2D 2Dmapppm/s). map ofoftheAinstantaneous thehand-drawn instantaneous ellipse outlines the concentration concentration approximate (ppm), (ppm), with withodor odor plume hits shape based on the location of bouts. thr = 0.04 The average bout frequency along the y-axis is shown in hitsrepresented representedbybyblue bluecircles circles(b(𝑏 0.04ppm/s). ppm/s).AAhand-drawn hand-drawnellipse ellipseoutlines outlinesthetheapproximate approximate the panel on plume the left. the The box plotsofbelow the map represents the instantaneous the concentration along plumeshape shapebased basedon on thelocation location ofbouts. bouts.The Theaverage averageboutboutfrequency frequencyalong along they-axis y-axisisisshown showninin the x-axis; (B) Drone trajectory in the z-axis. (C) Temporal evolution of the panel on the left. The box plots below the map represents the instantaneous concentration alongthe instantaneous the panel on the left. The box plots below the map represents the instantaneous concentration along the thex-axis; (B) (B) x-axis; Drone trajectory Drone in the z-axis. trajectory in the(C) Temporal z-axis. (C) evolution Temporalofevolution the instantaneous of the concentration instantaneous Sensors 2019, 19, 478 18 of 25 Sensors 2018, 18, x FOR PEER REVIEW 18 of 26 concentration (ppm), (ppm),bouts with detected with detected highlightedbouts highlighted in red (the blackinstar red indicates (the blackthe star indicates start the start of the bout). Theof bout the bout). The bout frequency (gray line) is computed using a sliding window of 5 s. The identifiers frequency (gray line) is computed using a sliding window of 5 s. The identifiers R1–R4 between panels R1– R4and (B) between panels the (C) indicate (B) area and of (C)the indicate map inthe areathe which of drone the map in which is flying the moment. at each drone is The flying at each maximum moment. The maximum instantaneous concentration and the maximum bout frequency instantaneous concentration and the maximum bout frequency are indicated by a green star and a blue are indicated by a green triangle, star and a blue triangle, respectively. respectively. Figure 19. Results of Experiment 2 when bthr is increased to 0.18 ppm/s. (top) 2D map of the Figure 19. Results instantaneous of Experiment concentration (ppm), 2with whenodor𝑏 hitsisrepresented increased toby0.18 blueppm/s. circles. (top) 2D map of A hand-drawn the ellipse instantaneous concentration (ppm), with odor hits represented by blue circles. A hand-drawn outlines the approximate plume shape based on the location of bouts. (bottom) Temporal evolution of ellipse outlines the the approximate instantaneous plume(ppm), concentration shapewith based on the location detected of bouts. (bottom) bouts highlighted in red (theTemporal black starevolution indicates the start of the bout). The bout frequency (gray line) is computed using a sliding (the of the instantaneous concentration (ppm), with detected bouts highlighted in red blackofstar window 5 s. indicates The maximumthe start boutoffrequency the bout). is The bout frequency indicated by a blue (gray line) is computed using a sliding window of triangle. 5 s. The maximum bout frequency is indicated by a blue triangle. 3.5. Experiment 3: Localization of a Source Hidden Inside a Power Outlet Box (h = 0.9 m) 3.5. Experiment 3: Localization of a Source Hidden Inside a Power Outlet Box (h = 0.9 m) In this experiment, the gas source was placed inside a power outlet box and a bubbler was used In thisthe to increase experiment, the release rate (nogas source was artificial placed airflow). The inside a power outlet gas distribution map box and a 20A) (Figure bubbler wasaused shows local to increase the accumulation ofrelease ratenear gas in R3 (no the artificial source airflow). The gas distribution (max concentration = 8 ppm), map (Figurein20A) whereas shows the rest a local of the map accumulation of gas in R3 near the source (max concentration = 8 ppm), whereas the gas concentration is below 1 ppm. The bout map does not resemble the elongated plume observed in the rest of the map the gas concentration is below 1 ppm. The bout map does not resemble in previous experiments probably due to the absence of induced airflow. During the first part of thethe elongated plume observed inatprevious exploration h = 0.9 experiments m, the measured probably due to thewas concentration absence of induced constantly below airflow. the LOD During and nothebouts first part of the exploration at ℎ 0.9 m, the measured concentration was constantly were detected (Figure 20C), which represents a clear challenge for a robot attempting to reactively below the LOD and no bouts were detected (Figure 20C), which represents a clear challenge for track the plume. Only when the drone approached the source (d < 3 m) the first bouts were detected a robot attempting to reactively (bouts #1–3).track The the plume.bouts remaining Only and when thethe drone approached maximum concentrationthewere sourcefound 3 m) the (𝑑 during the second first bouts pass were detected (bouts #1–3). The remaining bouts and the maximum concentration at h = 1.8 m, which may indicate a buoyant gas dispersion due to convective air currents produced were found during theheat by second insidepass ℎ 1.8 theatpower m, which outlet box. Themaylocation indicateofa the buoyant sourcegasin dispersion due to convective the instantaneous concentration air currents map produced and using by heat the bout inside thewere frequency power outlet very box. The accurate location in the of the(error x-y plane source ofin the 0.7 instantaneous m), but the bout concentration map and using the bout frequency were very accurate in the frequency yielded a lower localization error (0.77 m) than the instantaneous response (1.22 m) x-y plane (error of because 0.7 m), but the bout frequency yielded a lower localization error (0.77 m) than the instantaneous response the estimated height of the source in the first case (1.2 m) was closer to the true value (0.9 m) than in (1.22 m) because the estimated height of the source in the first case (1.2 m) was closer to the true value the second case (1.7 m). (0.9 m) than in the second case (1.7 m). Sensors 2019, 19, 478 19 of 25 Sensors 2018, 18, x FOR PEER REVIEW 19 of 26 Figure Figure 20. Results of 20. Results Experiment 3. of Experiment (A) 2D 3. (A) 2D map map ofof the the instantaneous instantaneous concentration concentration (ppm), (ppm), with with bouts bouts represented represented by byblue bluecircles circles (𝑏= 0.200.20 (bthr ppm/s). ppm/s).A hand-drawn A hand-drawn ellipse outlines ellipse the approximate outlines plume the approximate shape plumebased shape on the on based location of bouts. the location The average of bouts. bout bout The average frequency alongalong frequency the y-axis is shown the y-axis in the is shown in leftmost the leftmost panel. The box plot below the map represents the instantaneous concentration along the panel. The box plot below the map represents the instantaneous concentration along the x-axis; x-axis; (B) (B) Drone Drone trajectory trajectory in in the the z-axis. z-axis. (C) (C) Temporal Temporal evolution evolution of of the the instantaneous instantaneous concentration concentration (ppm), (ppm), with detected bouts highlighted in red (the black star indicates the start with detected bouts highlighted in red (the black star indicates the start of of the the bout). bout). The The bout bout frequency (gray line) is computed using a sliding window of 5 s. The identifiers R1–R4 frequency (gray line) is computed using a sliding window of 5 s. The identifiers R1-R4 between panels between panels (B) (B) and and (C) (C) indicate indicate the the area area of of the the map map inin which which the the drone drone is is flying flying atat each each moment. moment. The The maximum maximum of of the instantaneous concentration and the bout frequency are indicated by a green star and the instantaneous concentration and the bout frequency are indicated by a green star and aa blue blue triangle, triangle, respectively. respectively. 3.6. Overall Localization Results 3.6. Overall Localization Results A summary of the localization results is given in Table 2. The bout frequency with optimum noise A summary of the localization results is given in Table 2. The bout frequency with optimum threshold achieved comparable results to the instantaneous concentration, except in Experiment 2 noise threshold achieved comparable results to the instantaneous concentration, except in where the former achieved better performance (2.2 m versus 4.0 m error). The bout frequency ((µ + 3σ ) Experiment 2 where the former achieved better performance (2.2 m versus 4.0 m error). The bout threshold) produced relatively high localization errors. frequency (𝜇 + 3𝜎 threshold) produced relatively high localization errors. Table 2. Gas source localization error (m) in the three experiments, using the instantaneous concentration, the bout frequency with 𝜇 + 3𝜎 threshold or the bout frequency with optimum threshold. Experiment Instantaneous Bout Frequency Bout Frequency Sensors 2019, 19, 478 20 of 25 Table 2. Gas source localization error (m) in the three experiments, using the instantaneous concentration, the bout frequency with µ + 3σ threshold or the bout frequency with optimum threshold. Bout Frequency Experiment Instantaneous Concentration Bout Frequency (µ+3σ) (Optimum Threshold) 1 0.94 4.32 1.16 2 4.0 3.31 2.22 3 1.22 5.07 0.77 Mean 2.05 4.23 1.38 4. Discussion Our results suggest that a gas-sensitive NAV can be used for gas source localization and mapping in large indoor environments. The two GSL strategies compared in this paper require either building a map of the gas distribution or a map of the bouts. In contrast to previous works in which long measurement times are needed at predefined or adaptively chosen sampling locations, we have demonstrated that a rough approximation of such maps can be obtained in very short time with concentration measurements acquired in motion. Both maps seem coherent with respect to the true source location and wind direction, and not only enable the detection of the source with relatively small localization errors but also provide a rich visual interpretation of the gas distribution, especially if the bouts are overlaid on the gas distribution map. We adapted the bout detection algorithm for real-time operation and optimized the noise threshold, discussing how it affects the results in each experiment. The instantaneous gas distribution provides a more intuitive representation of the gas source location whereas the bout map can outline the gas plume without requiring anemometry. The experiments presented in this paper also demonstrate that the air flow generated by the propellers greatly affects the gas sensor signals, but key information for GSL can still be extracted by signal processing methods. This appears to contradict the results reported by Rossi et al. [29] who claimed that the signals of a MOX sensor on board of a Crazyflie 2.0 are useless if the propellers are switched on. It is difficult to reason about the causes of this discrepancy because Rossi et al. do not show nor describe the experimental arena in which their experiments were performed, including the type of gas source and its location relative to the drone position. If they attempted to measure the gas coming out from an evaporating chemical source placed right below the drone, we have seen in this work that the downwash of the drone prevents most of the gas reaching the sensor surface and this produces raw signals with much less fluctuations than they would have if the same measurements were carried out with the propellers turned off (see green trace in Figure 12). We have also seen that such characteristic fluctuations can be still extracted from the smoothed derivative of the response (Figure 13). Another possible reason for this difference may be that different MOX sensors were used (MiCS 5525 for CO detection in [29] and Figaro TGS 8100 in this work) or they were operated differently. While we operate the sensor in continuous mode, Rossi et al. apply duty-cycling to the sensor heater to reduce the power consumption, which is known to degrade the signal quality [54]. Other reasons may be the signal processing carried out, the test gas concentrations or the range of variation of environmental conditions. If the gas concentrations are close to the LOD of the sensor, any turbulence created by the propellers or slight changes in ambient temperature or humidity may hinder distinguishing the response to the gas from the noise [72]. We dedicated a substantial portion of the effort to calibrate the MOX sensor in similar conditions to the test scenario and to find a suitable linear calibration model that enables the computation of the LOD using the standard equations. The log-log transformation linearized the MOX response in the range 0–50 ppm and the LOD was estimated as 1.1 ppm. This is something not usually done in the field of mobile robot olfaction since calibration of gas sensors is time-consuming and requires expensive dedicated equipment. The sensor signals are often uncalibrated [73–76] or only scaled to unit range [77–80], which produces maps that represent the sensor response instead of the absolute concentration. Response maps are not easily interpretable due to the non-linear sensor response Sensors 2019, 19, 478 21 of 25 and are sensor-dependent; i.e., two response maps obtained with different gas sensors may show notable differences even if the underlying gas distribution is similar. Scaling is a linear operation that cannot remove non-linearities of the original signals. Some authors assume a linear response within a certain concentration range [75,76]; however, for most gas sensing technologies (e.g., MOX, gasFETs or thermoelectric sensors) this is only true within a narrow range of concentrations (e.g., 0–10 ppm) often exceeded by the typical concentrations obtained in real experiments. Most GSL algorithms require a detection threshold that is used to differentiate when a given sensor reading is produced by the gas plume or by background noise. Algorithms based on ‘odor hits’ require such a threshold to binarize the sensor response or to filter out low-amplitude bouts produced by noise. The rapid decay in chemical concentration with increasing distance from the source combined with incorrectly estimated detection thresholds may lead to a reduced area where the plume can effectively be detected (if the threshold is too high) or a high number of false alarms (if the threshold is too low). This has not been reported in the literature as an issue in GSL experiments because the concentration was often unrealistically high and the arena too small, i.e., the robot was very close to the source during the whole experiment. In this work, we challenged the drone by placing the source up to 17 m from the starting point or hiding it in the suspended ceiling or inside a power outlet box, leading to concentrations near the LOD of the sensor. It seems that for proper estimation of the noise threshold it is critical to have some prior information regarding the source intensity. A threshold computed exclusively from measurements in the absence of gas (e.g., µ + 3σ threshold) may only work when the concentrations in the room are small. If the background concentration in the room increases, for example due to gas dispersed by the propellers of the drone, the threshold must be increased to reduce false bout detections. This suggests a multi-threshold approach, i.e., building bout maps with different thresholds and searching for consistent source predictions over multiple threshold values. An obvious future direction for nano-drones is to achieve autonomous navigation, i.e., without resorting to external localization systems. This requires three-dimensional obstacle detection and self-localization, which is not an easy task whatsoever because there are many degrees of freedom and the low payload does often not allow for range sensors. In this line, the Fast Lightweight Autonomy (FLA) program (https://www.darpa.mil/program/fast-lightweight-autonomy) from the Defense Advanced Research Projects Agency (DARPA) is exploring novel perception methods and algorithms for small autonomous UAVs to fly at speeds up to 20 m/s in cluttered environments with no remote pilot, no communication links, no GPS guidance and no pre-programmed map of the area. In phase 2 of the program (June 2018), successful tests were performed with drones weighing around 2 kg equipped with a single camera, flying through a narrow window into a building, searching rooms, creating a 3-D map of the interior and exiting the building through an open doorway. In this paper we have demonstrated that a gas sensor could complement the camera in this application, adding a gas distribution layer to the acquired map of the environment which could be eventually analyzed by an operator to identify potential leak sources or victims. Before reaching this goal, the SNAV shall be tested in more complex scenarios, for example in environments with chemical interferences or multiple gas sources. In this context, it may be required to use multiple MOX sensors, integrate air flow information or implement more advanced localization/mapping algorithms, for example bioinspired reactive plume tracking or BASED maps of bouts [81,82] (i.e., using bouts as events or as weighted events). Author Contributions: Conceptualization, J.B. and S.M.; methodology, J.B. and V.H.; software, J.B.; validation, V.H., A.L. and S.M.; formal analysis, J.B.; investigation, J.B.; resources, A.L. and S.M.; data curation, J.B.; writing—original draft preparation, J.B.; writing—review and editing, all authors; visualization, J.B.; supervision, A.L. and S.M. Funding: This research was funded by Spanish MINECO, grant numbers BES-2015-071698 (Severo-Ochoa) and TEC2014-59229-R (SIGVOL). Sensors 2019, 19, 478 22 of 25 Acknowledgments: CERCA Programme/Generalitat de Catalunya. The Signal and Information Processing for Sensor Systems group is a consolidated Grup de Recerca de la Generalitat de Catalunya and has support from the Departament d’Universitats, Recerca i Societat de la Informació de la Generalitat de Catalunya (expedient 2014-SGR-1445). Conflicts of Interest: The authors declare no conflict of interest. References 1. Ma, K.Y.; Chirarattananon, P.; Fuller, S.B.; Wood, R.J. Controlled flight of a biologically inspired, insect-scale robot. Science 2013, 340, 603–607. [CrossRef] 2. Hassanalian, M.; Abdelkefi, A. Classifications, applications, and design challenges of drones: A review. Prog. Aerosp. Sci. 2017, 91, 99–131. [CrossRef] 3. Everaerts, J. The use of unmanned aerial vehicles (uavs) for remote sensing and mapping. Int. Arch. Photogramm. Remote Sens. Spat. Inf. Sci. 2008, 37, 1187–1192. 4. Khan, A.; Schaefer, D.; Tao, L.; Miller, D.J.; Sun, K.; Zondlo, M.A.; Harrison, W.A.; Roscoe, B.; Lary, D.J. Low power greenhouse gas sensors for unmanned aerial vehicles. Remote Sens. 2012, 4, 1355–1368. [CrossRef] 5. Berman, E.S.F.; Fladeland, M.; Liem, J.; Kolyer, R.; Gupta, M. Greenhouse gas analyzer for measurements of carbon dioxide, methane, and water vapor aboard an unmanned aerial vehicle. Sens. Actuators B Chem. 2012, 169, 128–135. [CrossRef] 6. Carrozzo, M.; De Vito, S.; Esposito, E.; Salvato, M.; Formisano, F.; Massera, E.; Di Francia, G.; Veneri, P.D.; Iadaresta, M.; Mennella, A. UAV Intelligent Chemical Multisensor Payload for Networked and Impromptu Gas Monitoring Tasks. In Proceedings of the 2018 5th IEEE International Workshop on Metrology for AeroSpace (MetroAeroSpace), Rome, Italy, 20–22 June 2018; pp. 112–116. 7. Chang, C.C.Y.; Wang, J.L.; Chang, C.C.Y.; Liang, M.C.; Lin, M.R. Development of a multicopter-carried whole air sampling apparatus and its applications in environmental studies. Chemosphere 2016, 144, 484–492. [CrossRef] 8. Xie, T.; Liu, R.; Hai, R.T.; Hu, Q.H.; Lu, Q. UAV platform based atmospheric environmental emergency monitoring system design. J. Appl. Sci. 2013, 13, 1289–1296. [CrossRef] 9. Rossi, M.; Brunelli, D. Autonomous Gas Detection and Mapping With Unmanned Aerial Vehicles. IEEE Trans. Instrum. Meas. 2016, 65, 765–775. [CrossRef] 10. McGonigle, A.J.S.; Aiuppa, A.; Giudice, G.; Tamburello, G.; Hodson, A.J.; Gurrieri, S. Unmanned aerial vehicle measurements of volcanic carbon dioxide fluxes. Geophys. Res. Lett. 2008, 35. [CrossRef] 11. Shinohara, H. Composition of volcanic gases emitted during repeating Vulcanian eruption stage of Shinmoedake, Kirishima volcano, Japan. Earth Planets Space 2013, 65, 667–675. [CrossRef] 12. Rüdiger, J.; Tirpitz, J.L.; de Moor, J.M.; Bobrowski, N.; Gutmann, A.; Liuzzo, M.; Ibarra, M.; Hoffmann, T. Implementation of electrochemical, optical and denuder-based sensors and sampling techniques on UAV for volcanic gas measurements: Examples from Masaya, Turrialba and Stromboli volcanoes. Atmos. Meas. Tech. 2018, 11, 2441–2457. [CrossRef] 13. Mori, T.; Hashimoto, T.; Terada, A.; Yoshimoto, M.; Kazahaya, R.; Shinohara, H.; Tanaka, R. Volcanic plume measurements using a UAV for the 2014 Mt. Ontake eruption the Phreatic Eruption of Mt. Ontake Volcano in 2014 5. Volcanology. Earth Planets Space 2016, 68. [CrossRef] 14. Astuti, G.; Giudice, G.; Longo, D.; Melita, C.D.; Muscato, G.; Orlando, A. An overview of the “volcan project”: An UAS for exploration of volcanic environments. J. Intell. Robot. Syst. Theory Appl. 2009. [CrossRef] 15. Neumann, P.P.; Kohlhoff, H.; Hüllmann, D.; Lilienthal, A.J.; Kluge, M. Bringing Mobile Robot Olfaction to the next dimension—UAV-based remote sensing of gas clouds and source localization. In Proceedings of the 2017 IEEE International Conference on Robotics and Automation (ICRA), Singapore, 29 May–3 June 2017; pp. 3910–3916. 16. Golston, L.M.; Aubut, N.F.; Frish, M.B.; Yang, S.; Talbot, R.W.; Gretencord, C.; McSpiritt, J.; Zondlo, M.A. Natural gas fugitive leak detection using an unmanned aerial vehicle: Localization and quantification of emission rate. Atmosphere 2018, 9, 333. [CrossRef] 17. Krüll, W.; Tobera, R.; Willms, I.; Essen, H.; Von Wahl, N. Early forest fire detection and verification using optical smoke, gas and microwave sensors. Procedia Eng. 2012, 45, 584–594. [CrossRef] Sensors 2019, 19, 478 23 of 25 18. Merino, L.; Caballero, F.; Martínez-de Dios, J.R.; Ferruz, J.; Ollero, A. A cooperative perception system for multiple UAVs: Application to automatic detection of forest fires. J. F. Robot. 2006, 23, 165–184. [CrossRef] 19. Pfeifer, J.; Khanna, R.; Constantin, D.; Popovic, M.; Galceran, E.; Walter, A.; Siegwart, R.; Liebisch, F. Towards automatic UAV data interpretation. In Proceedings of the International Conference of Agricultural Engineering 2016, At Aahus, Denmark, 26–29 June 2016. 20. Roldán, J.J.; Joossen, G.; Sanz, D.; del Cerro, J.; Barrientos, A. Mini-UAV based sensory system for measuring environmental variables in greenhouses. Sensors 2015, 15, 3334–3350. [CrossRef] 21. Pobkrut, T.; Eamsa-Ard, T.; Kerdcharoen, T. Sensor drone for aerial odor mapping for agriculture and security services. In Proceedings of the 2016 13th International Conference on Electrical Engineering/Electronics, Computer, Telecommunications and Information Technology (ECTI-CON), Chiang Mai, Thailand, 28 June–1 July 2016. [CrossRef] 22. Lega, M.; Napoli, R.M.A. A new approach to solid waste landfills aerial monitoring. WIT Trans. Ecol. Environ. 2008, 109, 193–199. [CrossRef] 23. Allen, G.; Hollingsworth, P.; Kabbabe, K.; Pitt, J.R.; Mead, M.I.; Illingworth, S.; Roberts, G.; Bourn, M.; Shallcross, D.E.; Percival, C.J. The development and trial of an unmanned aerial system for the measurement of methane flux from landfill and greenhouse gas emission hotspots. Waste Manag. 2018. [CrossRef] 24. Emran, B.J.; Tannant, D.D.; Najjaran, H. Low-altitude aerial methane concentration mapping. Remote Sens. 2017, 9, 823. [CrossRef] 25. Daniel, K.; Dusza, B.; Lewandowski, A.; Wietfeld, C. Airshield: A system-of-systems muav remote sensing architecture for disaster response. In Proceedings of the 2009 3rd Annual IEEE Systems Conference, Vancouver, BC, Canada, 23–26 March 2009; pp. 196–200. [CrossRef] 26. Alvarado, M.; Gonzalez, F.; Fletcher, A.; Doshi, A. Towards the development of a low cost airborne sensing system to monitor dust particles after blasting at open-pit mine sites. Sensors 2015, 15, 19667–19687. [CrossRef] [PubMed] 27. Villa, T.; Gonzalez, F.; Miljievic, B.; Ristovski, Z.; Morawska, L. An Overview of Small Unmanned Aerial Vehicles for Air Quality Measurements: Present Applications and Future Prospectives. Sensors 2016, 16, 1072. [CrossRef] [PubMed] 28. Pajares, G. Overview and Current Status of Remote Sensing Applications Based on Unmanned Aerial Vehicles (UAVs). Photogramm. Eng. Remote Sens. 2015, 81, 281–330. [CrossRef] 29. Rossi, M.; Brunelli, D. Gas Sensing on Unmanned Vehicles: Challenges and Opportunities. In Proceedings of the 2017 New Generation of CAS (NGCAS), Genova, Italy, 6–9 September 2017; pp. 117–120. 30. Fahad, H.M.; Shiraki, H.; Amani, M.; Zhang, C.; Hebbar, V.S.; Gao, W.; Ota, H.; Hettick, M.; Kiriya, D.; Chen, Y.; et al. Room temperature multiplexed gas sensing using chemical-sensitive 3. 5-nm-thin silicon transistors. Sci. Adv. 2017, 3, 1–9. [CrossRef] [PubMed] 31. Dunkley, O.; Engel, J.; Sturm, J.; Cremers, D. Visual-Inertial Navigation for a Camera-Equipped 25 g Nano-Quadrotor. In Proceedings of the IROS Aerial Open Source Robotics Workshop, Chicago, IL, USA, 14–18 September 2014; pp. 4–5. 32. Preiss, J.A.; Honig, W.; Sukhatme, G.S.; Ayanian, N. Crazyswarm: A large nano-quadcopter swarm. In Proceedings of the 2017 IEEE International Conference on Robotics and Automation (ICRA), Singapore, 29 May–3 June 2017; pp. 3299–3304. 33. Farid, Z.; Nordin, R.; Ismail, M. Recent advances in wireless indoor localization techniques and system. J. Comput. Networks Commun. 2013, 2013, 185138. [CrossRef] 34. Monroy, J.G.; González-Jiménez, J.; Blanco, J.L. Overcoming the slow recovery of MOX gas sensors through a system modeling approach. Sensors 2012, 12, 13664–13680. [CrossRef] 35. Lilienthal, A.; Duckett, T. Building gas concentration gridmaps with a mobile robot. Rob. Auton. Syst. 2004, 48, 3–16. [CrossRef] 36. Kowadlo, G.; Russell, R.A. Robot odor localization: A taxonomy and survey. Int. J. Rob. Res. 2008, 27, 869–894. [CrossRef] 37. Lochmatter, T. Bio-Inspired and Probabilistic Algorithms for Distributed Odor Source Localization using Mobile Robots; École polytechnique fédérale de Lausanne (EPFL): Lausanne, Switzerland, 2010; Volume 4628. 38. Hernandez Bennetts, V.; Lilienthal, A.J.; Neumann, P.P.; Trincavelli, M. Mobile Robots for Localizing Gas Emission Sources on Landfill Sites: Is Bio-Inspiration the Way to Go? Front. Neuroeng. 2012, 4. [CrossRef] Sensors 2019, 19, 478 24 of 25 39. Lilienthal, A.; Duckett, T. Experimental analysis of gas-sensitive Braitenberg vehicles. Adv. Robot. 2004, 18, 817–834. [CrossRef] 40. Ishida, H.; Nakamoto, T.; Moriizumi, T. Remote sensing of gas/odor source location and concentration distribution using mobile system. Sens. Actuators B Chem. 1998, 49, 52–57. [CrossRef] 41. Pang, S.; Farrell, J.A. Chemical plume source localization. IEEE Trans. Syst. Man Cybern. Part B Cybern. 2006, 36, 1068–1080. [CrossRef] 42. Vergassola, M.; Villermaux, E.; Shraiman, B.I. “Infotaxis” as a strategy for searching without gradients. Nature 2007, 445, 406–409. [CrossRef] [PubMed] 43. Pomareda, V.; Magrans, R.; Jiménez-Soto, J.M.; Martínez, D.; Tresánchez, M.; Burgués, J.; Palacín, J.; Marco, S. Chemical source localization fusing concentration information in the presence of chemical background noise. Sensors 2017, 17, 904. [CrossRef] [PubMed] 44. Turner, D.B. Workbook of Atmospheric Dispersion Estimates: An Introduction to Dispersion Modeling; CRC Press: Boca Raton, FL, USA, 1994. 45. Farrell, J.A.; Murlis, J.; Long, X.; Li, W.; Cardé, R.T. Filament-based atmospheric dispersion model to achieve short time-scale structure of odor plumes. Environ. Fluid Mech. 2002, 2, 143–169. [CrossRef] 46. Thomson, D.J. Criteria for the selection of stochastic models of particle trajectories in turbulent flows. J. Fluid Mech. 1987, 180, 529–556. [CrossRef] 47. Falkovich, G.; Gawȩdzki, K.; Vergassola, M. Particles and fields in fluid turbulence. Rev. Mod. Phys. 2001, 73, 913–975. [CrossRef] 48. Sutton, O.G. The problem of diffusion in the lower atmosphere. Q. J. R. Meteorol. Soc. 1947, 73, 257–281. [CrossRef] 49. Bakkum, E.A.; Duijm, N.J. Vapour Cloud Dispersion; CPR E: London, UK, 1997; Volume 14. 50. Luo, B.; Meng, Q.H.; Wang, J.Y.; Sun, B.; Wang, Y. Three-dimensional gas distribution mapping with a micro-drone. In Proceedings of the 2015 34th Chinese Control Conference (CCC), Hangzhou, China, 28–30 July 2015; pp. 6011–6015. [CrossRef] 51. Lilienthal, A.J.; Reggente, M.; Trinca, M.; Blanco, J.L.; Gonzalez, J. A statistical approach to gas distribution modelling with mobile robots—The Kernel DM+V algorithm. In Proceedings of the 2009 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), St. Louis, MO, USA, 10–15 October 2009; pp. 570–576. 52. Hayes, A.T.; Martinoli, A.; Goodman, R.M. Distributed odor source localization. IEEE Sens. J. 2002, 2, 260–271. [CrossRef] 53. Neumann, P.P.; Asadi, S.; Lilienthal, A.J.; Bartholmai, M.; Schiller, J.H. Autonomous gas-sensitive microdrone: Wind vector estimation and gas distribution mapping. IEEE Robot. Autom. Mag. 2012, 19, 50–61. [CrossRef] 54. Burgués, J.; Marco, S. Low power operation of temperature-modulated metal oxide semiconductor gas sensors. Sensors 2018, 18, 339. [CrossRef] [PubMed] 55. Lilienthal, A.; Zell, A.; Wandel, M.; Weimar, U. Sensing odour sources in indoor environments without a constant airflow by a mobile robot. In Proceedings of the 2001 IEEE International Conference on Robotics and Automation (ICRA), Seoul, Korea, 21–26 May 2001; Volume 4, pp. 1–6. 56. Atema, J. Chemical signals in the marine environment: Dispersal, detection, and temporal signal analysis. Proc. Natl. Acad. Sci. USA 1995, 92, 62–66. [CrossRef] [PubMed] 57. Farah, A.; Duckett, T. Reactive Localisation of an Odour Source by a learning Mobile Robot. In Proceedings of the Second Swedish Workshop on Autonomous Robotics, Stockholm, Sweden, 11–12 October 2002; pp. 29–38. 58. Weissburg, M.J.; Dusenbery, D.B.; Ishida, H.; Janata, J.; Keller, T.; Roberts, P.J.W.; Webster, D.R. A multidisciplinary study of spatial and temporal scales containing information in turbulent chemical plume tracking. Environ. Fluid Mech. 2002, 2, 65–94. [CrossRef] 59. Webster, D.R.; Weissburg, M.J. Chemosensory guidance cues in a turbulent chemical odor plume. Limnol. Oceanogr. 2001, 46, 1034–1047. [CrossRef] 60. Gonzalez-Jimenez, J.; Monroy, J.G.; Blanco, J.L. The multi-chamber electronic nose-an improved olfaction sensor for mobile robotics. Sensors 2011, 11, 6145–6164. [CrossRef] [PubMed] 61. Batog, P.; Wołczowski, A. Odor markers detection system for mobile robot navigation. Procedia Eng. 2012, 47, 1442–1445. [CrossRef] Sensors 2019, 19, 478 25 of 25 62. Marco, S.; Pardo, A.; Davide, F.A.M.; Di Natale, C.; D’Amico, A.; Hierlemann, A.; Mitrovics, J.; Schweizer, M.; Weimar, U.; Göpel, W. Different strategies for the identification of gas sensing systems. Sens. Actuators B Chem. 1996, 34, 213–223. [CrossRef] 63. Pardo, A.; Marco, S.; Samitier, J.; Davide, F.A.M.; Di Natale, C.; D’Amico, A. Dynamic measurements with chemical sensor arrays based on inverse modelling. In Proceedings of the IEEE Instrumentation and Measurement Technology Conference, Brussels, Belgium, 4–6 June 1996; Volume 1. 64. Fonollosa, J.; Sheik, S.; Huerta, R.; Marco, S. Reservoir computing compensates slow response of chemosensor arrays exposed to fast varying gas concentrations in continuous monitoring. Sens. Actuators B Chem. 2015, 215, 618–629. 65. Schmuker, M.; Bahr, V.; Huerta, R. Exploiting plume structure to decode gas source distance using metal-oxide gas sensors. Sens. Actuators B Chem. 2016, 235, 636–646. [CrossRef] 66. Bitcraze, A.B. Getting Started with the Loco Positioning System. Available online: https://www.bitcraze.io/ getting-started-with-the-loco-positioning-system/ (accessed on 7 July 2018). 67. DecaWave DWM1000 Datasheet; DecaWave: Dublin, Ireland, 2016. 68. Nelson, G. Gas Mixtures: Preparation and Control; CRC Press: Boca Raton, FL, USA, 1992. 69. Currie, L.A. Nomenclature in evaluation of analytical methods including detection and quantification capabilities (IUPAC Recommendations 1995). Pure Appl. Chem. 1995, 67, 1699–1723. [CrossRef] 70. Burgués, J.; Jimenez-Soto, J.M.; Marco, S. Estimation of the limit of detection in semiconductor gas sensors through linearized calibration models. Anal. Chim. Acta 2018, 1013, 13–25. [CrossRef] [PubMed] 71. Pashami, S.; Lilienthal, A.J.; Trincavelli, M. Detecting changes of a distant gas source with an array of MOX gas sensors. Sensors 2012, 12, 16404–16419. [CrossRef] [PubMed] 72. Burgués, J.; Marco, S. Multivariate estimation of the limit of detection by orthogonal partial least squares in temperature-modulated MOX sensors. Anal. Chim. Acta 2018, 1019, 49–64. [CrossRef] [PubMed] 73. Shakaff, A.; Yeon, A.; Kamarudin, K.; Visvanathan, R. Gas Source Localization via Behaviour Based Mobile Robot and Weighted Arithmetic Mean Gas Source Localization via Behaviour Based Mobile Robot and Weighted Arithmetic Mean. IOP Conf. Ser. Mater. Sci. Eng. 2018, 318, 012049. [CrossRef] 74. Li, J.G.; Sun, B.; Zeng, F.L.; Liu, J.; Yang, J.; Yang, L. Experimental study on multiple odor sources mapping by a mobile robot in time-varying airflow environment. In Proceedings of the Chinese Control Conference (CCC), Chengdu, China, 27–29 July 2016; pp. 6032–6037. 75. Kowadlo, G.; Russell, R.A. Using naïve physics for odor localization in a cluttered indoor environment. Auton. Robots 2006, 20, 215–230. [CrossRef] 76. Lilienthal, A.J.; Loutfi, A.; Duckett, T. Airborne chemical sensing with mobile robots. Sensors 2006, 6, 1616–1678. [CrossRef] 77. Reggente, M.; Lilienthal, A.J. Three-dimensional statistical gas distribution mapping in an uncontrolled indoor environment. AIP Conf. Proc. 2009, 1137, 109–112. 78. Pashami, S.; Lilienthal, A.J.; Schaffernicht, E.; Trincavelli, M. TREFEX: Trend estimation and change detection in the response of MOX gas sensors. Sensors 2013, 13, 7323–7344. [CrossRef] 79. Loutfi, A.; Coradeschi, S.; Lilienthal, A.J.; Gonzalez, J. Gas distribution mapping of multiple odour sources using a mobile robot. Robotica 2009, 27, 311–319. [CrossRef] 80. Lilienthal, A.; Reimann, D.; Zell, A. Gas Source Tracing with a Mobile Robot Using an Adapted Moth Strategy. Auton. Mob. Syst. 2003, 150–160. [CrossRef] 81. Lilienthal, A.; Trincavelli, M.; Schaffernicht, E. It’s always smelly around here! Modeling the Spatial Distribution of Gas Detection Events with BASED Grid Maps. In Proceedings of the 15th International Symposium on Olfaction and Electronic Nose (ISOEN 2013), Daegu, Korea, 2–5 July 2013; Volume 27. 82. Schaffernicht, E.; Trincavelli, M.; Lilienthal, A.J. Bayesian Spatial Event Distribution Grid Maps for Modeling the Spatial Distribution of Gas Detection Events. Sens. Lett. 2014, 12, 1142–1146. [CrossRef] © 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).