Gamma ray measurements involved in monitoring technologies of field conditions are of vital importance for environmental safety and radiation protection. This chapter addresses the method of cooperative gamma sensing using multiple unmanned aerial systems. Section 1 provides an introduction. The design of semiconductor and scintillation gamma ray sensors integrated into aerial robotic platforms is discussed in Section 2, along with the fusion of time-stamped radiation data with position information using the real-time kinematic positioning technique. Section 3 addresses the multirobot contour mapping of radiation fields. Computer simulation of radiation contour mapping is discussed in Section 4. Experimental verification of the contour mapping and source-seeking algorithm is described in Section 5. Section 6 summarizes results of the project.
- gamma rays
- remote sensing
- unmanned aerial system
- robot operating system
- contour mapping
Technologies of gamma ray sensing in radiation zones are vital for environmental safety and radiation safety. Hazardous radiological wastes can be a result of large nuclear projects; radiological materials can be displaced, lost, or smuggled [1, 2, 3, 4]. Additionally, they can be released into the environment due to man-made accidents or natural disasters. One example is the Fukushima Daiichi nuclear power plant disaster  that led to radiological contamination of the plant’s infrastructure and territory adjacent to the plant. The use of Unmanned Aerial Systems (UAS) allows for the remote radiation sensing, mapping, and wide area search operations while keeping users away from the risk of exposure.
Small scale, multiuse UAS platforms equipped with navigation and radiation sensing capabilities permit for the surveillance in hard to reach, hazardous areas while allowing for the measurements to be dynamically tracked and mapped . The measured data could be used for situational awareness and further analysis of radiation fields in temporal and space domains. Furthermore, based on the cooperative sensing algorithms, the UAS swarm can be programmed to search for unattended radiation sources. In order to accomplish radiation monitoring tasks using UAS, gamma ray sensors should be integrated into a robotic aerial platform. The sensor’s data analysis should be automated and carried out onboard taking into account limited computational resources of flight computers.
2. Radiation sensors
Since UAS has weight and size limitations for onboard components so that the flight efficiency is not reduced due to the payload, the choice of radiation sensors is important. Radiation sensors should be swappable (attached and detached from the robotic platform) in the field conditions. The sensors should have low power requirements so that the integration into the platform does not affect the battery power. The data processing limitations should be also taken into account. To address this necessity, two types of flight ready ambient temperature radiation sensors were developed: a Cadmium Zinc Telluride (CZT) semiconductor sensor for high-resolution gamma spectrometry and the elpasolite scintillation sensor Cs2LiYCl6:Ce3+ (CLYC) for neutron and gamma measurements. Both sensors were designed as plug-and-play interchangeable modules.
2.1 CZT sensor
A CZT is a wide bandgap semiconductor . Gamma rays that interact with this detection medium deposit their energy into the CZT operating in a direct-conversion mode at an ambient temperature. Therefore, no cooling is needed which significantly simplifies the sensor’s integration into the robots and their practical use. These semiconductors are able to process more than 10 million photons per square millimeter per second.
The CZT module GR1-A (Kromek)  was integrated into a multicopter aerial platform. The GR1-A includes a 1 cm3 cubic CZT crystal. Electric signals generated by the CZT proportionally to energy of the incident gamma ray are amplified by a preamplifier and then processed by a shaping amplifier. A 4096 channel analyzer produces the discrete data array of the gamma spectrum that can be processed further. The scheme of the CZT sensor operation is shown in Figure 1 . The USB interface is used for the data communication as well as power. The power consumption of this module is about 250 milliwatts. The module’s volume is 2.5 × 2.5 × 6.1 cm3 and the weight is about 50 grams. The sensor’s energy range is from 30 keV to 3 MeV with the Full Width at Half Maximum of the peak (FWHM) energy resolution of less than 2% at photon energy of 662 keV. The module’s temperature range is from 0 to 40°C. The CZT sensor mounted onto the UAS platform (the DJI S1000+ octocopter) is shown in Figure 2 .
Robot Operating System (ROS)  was used for fusion of the data sent from multiple sensors connected to the robotic platform (radiation sensors and a GPS sensor). ROS is an open source tool, consisting of libraries used for robot applications. It allows including a number of independent nodes which communicate to each other, using a publishing and subscribing messages to the topics. Once a data array was sent from the CZT sensor via USB to an onboard Linux minicomputer (the Odroid model) and processed using ROS, a gamma ray spectrum is generated. The example of a measured spectrum using a 60Co gamma-ray source is shown in Figure 1 .
The gamma ray spectrum is then analyzed using a code to extract the peaks of interest along with their intensities. This allows for identification of isotopes (sources) that emit these gamma rays, and also for estimation of source strength. Computational algorithms are widely used for photon spectral analytics for various applications [10, 11, 12, 13, 14]. For the rapid spectral analysis on board of robotic platforms with a limited computing power, a simple, but robust algorithm is needed. The algorithm based on Mariscotti’s technique  was developed for the peak identification in the presence of a background in the measured spectrum. In this technique, peaks are approximated using a Gaussian function.
Within a small interval of several peak’s widths, which is applicable for CZT detectors and scintillation detectors, the spectrum’s background is represented as a linear polynomial. In this interval, number of counts in the spectrum channel is reconstructed as , where is a Gaussian function, and constants and are associated with the background. This method locates a peak where the 2nd derivative of the function is not zero.
The function includes the height of the peak that is centered at a channel . The value of is found as . For a discrete data , the 2nd derivative of was approximated as a centered finite-divided difference: , where is the channel’s width.
Since carries statistical errors, the data will also exhibit statistical fluctuations about expected values at each channel. Hence, a peak cannot be resolved in case of an expected value of which is equal to the standard deviation. For it corresponds to the minimum resolvable value of 600. To resolve much smaller peaks, function should be ‘smoothed’, consequently reducing the standard deviation. The smoothing was done by calculating an average of values over adjacent channels near the channel . This smoothing can be repeated times yielding , where are weighting factors. Its standard deviation is . Values and were determined for the CZT sensor based on the desired value.
This algorithm makes it possible to distinguish photopeaks from Compton shoulders in the spectral data. The peak centroid was found at the center of the Gaussian function. The intensity of the peak was calculated using the area under the Gaussian. The C code implementing this algorithm was written as a function within ROS. The spectrum data array can be processed using this function returning the energies of found peaks and their intensities, as a result minimizing the data transfer from the UAS platform to a ground control station.
2.2 CLYC sensor
A 2.54-cm diameter, 2.54-cm long cylindrical CLYC [16, 17] crystal was utilized in the design of a gamma-neutron sensor. The CLYC’s density is 3.31 g/cm3. Its scintillation light wavelength range is from 275 to 450 nm (the peak is at 370 nm). The crystal’s refractive index is 1.81 at 405 nm. This crystal was coupled with a super bialkali photomultiplier tube matching the CLYC’s emission wavelength range, a miniature digitizer, and a high voltage generator; all packaged in a custom housing ( Figure 3a ). The 6 Li isotope enrichment of the CLYC was 95%. Neutron detection was achieved via 6 Li(n,α)t reaction. The scintillation light yield of CLYC is 20,000 photons per 1 MeV for gamma rays and 70,000 photons per thermal neutron. CLYC’s scintillation properties allow for gamma spectrometry. This sensor operates without cooling. The measured gamma-ray FWHM energy resolution was 5% at 662 keV and 3.3% at 1332 keV.
Moreover, CLYC exhibits neutron-photon pulse shape discrimination (PSD) properties. The photon-induced excitation in the CLYC medium leads to a fast core to valence (CVL) decay and prompt cerium decay with 1 and 50 ns decay constants, respectively. A neutron event in CLYC causes a slow cerium self-trapped excitation (Ce-STE, 1000 ns decay constant) .
The digitized neutron and photon signals of the CLYC sensor are shown in Figure 3b . Each signal was analyzed using the eMorpho digitizer, generating three values saved in a list mode: the time stamp; the integral under the front of the signal’s curve assessed using the partial integration time; and the integral under the entire curve that is proportional to the energy of measured radiation. In order to segregate neutron signals from gamma-ray signals, a PSD value (calculated as ratio of the area under the tail part of the signal to the front part of the signal) was used. Because neutron signals have longer tails, it leads to larger PSD values than in the case of gamma-ray signals. The experimental plot of neutron-gamma PSD for the CLYC sensor using a PuBe source is shown in Figure 3c . Neutron and gamma ray signals are well separated in this plot. The neutron signals appear at the value of 3200 keV, electron equivalents (keVee). The figure of merit of neutron-gamma PSD for CLYC was evaluated as 2.3 using an approach described in .
2.3 Sensor integration with UAS
The ‘plug-and-play’ concept was used to integrate the radiation sensors into the UAS using ROS. This approach supports ‘hot swapping’ of the sensors into the UAS platform, meaning that the user does not need to set up sensor’s parameters each time the UAS is powered on. When the CZT or CLYC sensor is plugged, the operating system recognizes it, and then installs the associated driver to read and process the sensor’s data. The data is published and the message appears. Similarly, when the sensor is unplugged, the operating system terminates the process and deletes the sensor’s driver. The scheme of the plug-and-play operations with CZT and CLYC sensors is shown in Figure 4 .
The ROS programming environment was also utilized for fusion of the radiation sensor’s data with time and position data to enable the spatiotemporal analytics of measured radiation fields. To determine the UAS coordinates at a time of the gamma-ray measurement, the real-time kinematic (RTK) positioning technique was used. This navigation method enhances precision of the position data obtained using satellite positioning systems. Based on the measurements of the phase of the signal’s carrier wave, it employs a single reference ‘base’ station for real-time corrections of the UAS position. The scheme of the RTK GPS technique is shown in Figure 5 .
A base station contains a fixed Swift Duro RTK GPS receiver, having known coordinates. Multiple UAS platforms carry the RTK GPS and L1/L2 GPS antenna. Correction data are transmitted from the base station to the UAS’s RTK Piksi Multi GPS receiver thus allowing calculating the ionospheric error. This method permits raw data measurements with a frequency of 20 Hz. Accuracy of the UAS position estimation using the RTK technique is 10 mm horizontally and 15 mm vertically. It was experimentally verified using a single base station and four UAS platforms. Such precision of position estimation for multiple UAS platforms enables to use the RTK GPS in cooperative multirobot radiation monitoring tasks including contour mapping and source search.
3. Multirobot contour mapping of radiation fields
One of the approaches of cooperative multirobot sensing is the method of contour mapping of radiation areas. It is based on the use of several UAS platforms (the ‘swarm’) equipped with radiation sensors. This approach allows for the automatic determination of the contour in space that corresponds to a preset radiation dose. Thus, the multirobot system could locate and follow a boundary of the hazardous zone.
In this section, the contour mapping algorithm is presented along with a gradient direction estimation and heading angle calculation scheme for the swarm consisting of three UAS that are positioned in a circular formation in the two dimensional space. It is assumed that a gamma-ray sensor, CZT or CLYC, is mounted on each UAS platform. The gamma-ray data are time stamped and merged with the position data as discussed in the previous section.
3.1 Gradient direction estimation
The contour mapping is based on two components: the gradient direction estimation and the average radiation level calculated using the radiation measurement data from sensors mounted on the UAS platforms of the swarm. The average of a scalar field is estimated over a circular area of radius centered at a point as shown in Figure 6 . is the average radiation level calculated using the data from sensors of UASs flying in a circular formation:
Here, is the number of UAS platforms, is the intensity of the measured gamma peak of interest at a point by ith UAS. The formation center moves toward the increasing (a source-seeking method) or the constant (a contour mapping method) value of the average of sensor readings. To find the required direction of motion of the swarm’s center, the gradient of should be determined using multiple readings from the UAS’s sensors. It is assumed that measurements are taken out of the readings distributed inside the circle according to a known distribution (e.g., uniform or based on a model). Using the composite trapezoidal rule, the gradient is estimated [20, 21] as:
Here, and . It should be noted that the origin is moved to the center of the circle, and the integral is approximated by a finite number of measurements.
Three UAS platforms of the swarm () are equally distributed around a circle circumference. Horizontal and vertical components of the gradient can be calculated using Eq. (2). Then the formation center can be moved relative to its current position using the estimated direction of gradient for the source-seeking behavior for the swarm. In this technique, any number of UAS platforms can be used in the swarm. For accurate estimation of the gradient, these drones should be evenly distributed around a circle. In the gradient estimation algorithm, there is an inherent error in a scalar field. A relatively small change in distance can have a large effect on gamma measurement for each UAS depending on relative orientation of the swarm to the source as well as its distance from the source. In the contour mapping technique, three UAS platforms rotate about the swarm’s center to improve the gradient estimation by changing the relative direction of the radiation source with respect to the UAS. Figure 7a illustrates how a gradient estimation error is defined. It should be noted that increases as a distance to the source relative to the swarm’s center decreases. Therefore, the estimation error decreases when the source is located far away from the swarm ( Figure 7b ). The number of UAS in the swarm affects the gradient error as well. Figure 7c shows that the gradient estimation improves with more UAS in the swarm.
3.2 Heading angle
The bulk heading angle of the UAS swarm ψ depends on how far this swarm is from the desired radiation contour to be mapped. As the swarm approaches the desired contour, the heading angle must be directed to a tangential direction for this contour. When the swarm is far away from the contour, the heading angle will be directed toward the source, which is the source-seeking behavior of the swarm as shown in Figure 8 . Here, is an estimated steepest gradient direction, and ϕ is a control angle determining a bulk heading angle Ψ . The angles are measured with respect to a positive -axis. The control angle ϕ is determined by how far the desired contour is located based on the average radiation measurement: , where is the desired radiation intensity of the contour to be mapped and is an average radiation intensity measured using three UAS platforms. Here, an arbitrary constant of small magnitude is used along with to calculate a heading angle. Value of is determined based on the PID control action from the measurement difference with respect to the reference contour value :
where , , and are the proportional, integral, and derivative gains, respectively. As shown in Eq. (3), a heading angle ψ becomes 90 degrees for a large value of when the swarm is far away from the reference contour.
For the swarm that is near the reference contour, becomes small and ψ is close to zero. This leads to the equation for determining the heading angle:
Based on Eq. (4), the swarm will move in a tangential direction near the reference contour, and demonstrate a source-seeking behavior when it is far away from the source or the reference contour.
4. Computer simulation of radiation contour mapping
4.1 UAS swarm simulation in 1/R 2 field
To validate the contour mapping algorithm, two types of radiation fields were used. The first is based on a model for a radioactive source located at the distance from a sensor. This simplified model was employed to properly tune and improve contour mapping and source-seeking behaviors of the swarm. The algorithm was improved by including the UAS dynamics and adaptive spinning of the swarm (for the reduced flight trajectory of each UAS). The improved algorithm was validated using a realistic gamma field formed by multiple sources placed in the area with physical obstacles that was computed using the Monte Carlo Neutron and Particles code (MCNP) . MCNP is widely used for simulation of coupled photon, neutron, electron and particle transport in complex geometries [23, 24].
A stochastic nature of radiation sensing was taken into account in the sensor’s data. A random noise was introduced to the radius value of the model in the simulation: a noise of ±2.5 m was added, causing erratic gradient estimation and heading generation, as expected. To reduce this effect, a moving average filter was applied to both the gradient estimation and the heading angle generation. The simulation scheme is shown in Figure 9 , including the UAS spinning formation used to estimate the steepest gradient. The UAS dynamics was incorporated for the realistic simulation of flight trajectories of aerial platforms during contour mapping and source seeking. Instead of using just kinematic motions of each UAS for this simulation, their dynamics was used to calculate their flight trajectories. The effect of inclusion of dynamics of UAS platforms on their trajectories that the swarm is able to follow is shown in Figure 10 . A double integrator was incorporated in MATLAB to approximate the UAS dynamics. This provides more realistic simulation comparing to a case when only kinematics is considered.
A spinning formation of the swarm improves the estimation of radiation gradient direction [21, 25]. Spinning occurs around a virtual center of the UAS formation while the center of the formation moves along the desired vector to accomplish the contour mapping or source seeking. Figure 11a illustrates the swarm of three UAS platforms spinning around a virtual center to counteract error from the gradient estimation algorithm. A bold line shows a path traveled by the swarm’s center in mapping a contour without spinning ( Figure 11b ) and with spinning ( Figure 11c ). The non-spinning formation exhibits poor mapping performance primarily due to large gradient estimation errors which depend on the relative direction to the source with respect to the swarm.
Having the swarm circling around a virtual center, while this center travels in the direction dictated by the swarm heading algorithm, increases the total flight paths significantly. This led to development of an adaptive spin rate adjustment scheme to avoid unnecessary spinning formation when it is not required.
When the UAS swarm is far away from the source or the contour to be mapped, the direction to them is nearly fixed, thus it is not necessary to spin the formation. When the swarm is near the contour, it is necessary to spin the formation since most contours have a curvilinear shape. The criterion for spinning is chosen as a radius of curvature of the formation center’s path.
As shown in Figure 12a , the swarm starts to spin when it approaches the contour to be followed. In this simulation, a spin rate control was bound with a lower value of 0.05 rad/s and an upper value of 0.3 rad/s. There is a big spike in a radius of curvature plot due to the nearly linear path the swarm has to follow. During the source-seeking stage, a minimum spinning was maintained to save on the total flight paths. Total path length in terms of spin rates is shown in Figure 12b which demonstrates a nearly negligible difference in path length if the swarm spins clockwise or counterclockwise. Rotation rates of 0, ±0.075, ±0.1, ±0.2, and ± 0.3 rad/s were used to compare the average of the flight path lengths.
Multiple radiation sources were also considered. Figure 13 shows contour mapping simulation for three sources; the algorithm was capable to trace the desired contour reasonably well. Simulation was also done for a moving radiation source. The algorithm can accomplish tracing of the desired radiation contour well if the source moves reasonably slow. Figure 14 shows a moving source traveling at 0.07 meters per second along the line from a point at (10, 40) m to a point (40, 10) m. Mapping this source was possible in this particular case because the speed of the swarm was roughly seven times the speed of the source.
4.2 UAS swarm simulation in MCNP computed radiation field
A realistic gamma-ray flux distributions in a 3D volume with dimensions of 100 m × 100 m × 32 m that contained five photon sources with energy ranging from 3 to 6 MeV were computed using MCNP. A concrete building was used in the model to introduce a structure that attenuates and blocks radiation flux in certain areas of the volume monitored by onboard sensors of UAS. Therefore, complex gamma ray flux distributions were generated ( Figure 15 ).
The reference contour along with performance of the contour mapping algorithm overlaid onto this radiation field is shown in Figure 16 . A two dimensional radiation distribution at the height of 15 m was used in this simulation (the UAS platforms were moving in the same plane). In this simulation, the swarm’s starting position was located inside the radiation field at a point with coordinates (35 m; 35 m). The algorithm allowed mapping the desired contour with a reasonable accuracy.
Key algorithms of contour mapping and source seeking were validated in the indoor flight testbed. This testbed was outfitted with the OptiTrack motion capture technology which allows real-time feedback of the position and orientation of moving robots within the flight volume at 120 Hz rate.
The Crazyflie 2.0 by Birtcraze and the DJI Flamewheel 450 are two UAS platforms that were used; Crazyflie is an open source, easily modifiable, light weight quadcopter. Its small size enables validation of the contour mapping algorithm with a virtual source. The DJI Flamewheel 450 was used for validation of the source-seeking algorithm using a light source of the type. Each Flamewheel platform was equipped with an onboard omnidirectional light sensor. The Crazyflie and Flamewheel platforms with a single board computer mounted under the frame are shown in Figure 17a and Figure 17b , respectively.
5.1 Gradient estimation
To verify the algorithm of gradient estimation by three UAS in a circular formation, each platform was placed upside down roughly equiangular around the center of the flight volume as shown in Figure 18a . Each UAS was placed at 0.5 m radius mark away from the center. The light source, acting as a radiation source analog, was moved around in a circle concentric with the swarm. The data stream captured through a wireless sensor network was input to the gradient estimation algorithm. Positions of a light source and three UAS were identified by the OptiTrack motion capture system, and used to compute true values of the gradient based on the assumption.
Figure 18b shows that the gradient direction estimation scheme agrees with the measured data reasonably well. It should be noted that the source distance of 0.7 m was used, which is rather too small for actual contour mapping application in the field. As expected, gradient direction estimation error reached almost 30 degrees which is rather too big for accurate mapping operation. However, it should be noted that experimental measurement data have a good agreement with the computed data.
5.2 Source-seeking behavior
Due to space constraints of the indoor flight volume, the source-seeking experiment was carried out using multiple Flamewheel platforms to test the gradient estimation algorithm for the heading angle rather than the contour mapping. A light source was placed on a movable dolly within the flight volume, and moved along the -axis. It was shown that the gradient estimation algorithm was effective in determining the direction to the source using the light sensor’s data from each UAS.
The swarm was restricted to move along the -axis and the reference position generation was bound to ±0.5 m, in order to minimize risk of the swarm crashing into the walls of the flight volume. The tracking of the light source and three UAS platforms was carried out using the OptiTrack. Note that an embedded window to display source-seeking performance in real time using the motion capture system in the flight volume is shown in Figure 19a . Figure 19b shows how the UAS swarm’s center moves as the source is moved. As expected, the oscillatory motion of the swarm’s center occurs due to gradient estimation errors using a finite number of light sensors.
5.3 Contour mapping behavior
To demonstrate the effectiveness of the contour mapping algorithm in the indoor flight volume, three Crazyflie platforms were used. A virtual source was used due to a limited payload and communication capability of the Crazyflie UAS. As shown in Figure 20a , a virtual source was located on the ground and the OptitTrack tracked the source and each UAS in the swarm. The ‘source strength’ needed for the gradient estimation algorithm was calculated using the model where was obtained from the virtual source position data from each UAS. Figure 20b shows the experimental results on a plot of motion of the swarm’s center following the reference contour defined with respect to the virtual source located on the floor. It maps the reference contour within ±0.1 m, which is less than 8% of the size of the contour.
Ambient temperature CZT and CLYC sensors were integrated onto the UAS platform using the plug-and-play approach. The CZT sensor was designed for high resolution gamma spectroscopy. The CLYC sensor enables gamma and neutron measurements with an excellent neutron-gamma pulse shape discrimination with the figure of merit 2.3. The automated spectral analysis code locating peaks and calculating their intensities was developed for both sensors.
USB hardware connections were used to link the sensors and the main controller with the UAS power source. ROS was used for the data communication and data fusion. To streamline the process of bridging disparate components into a cohesive network, the collection of libraries describing the publisher/subscriber communication of ROS nodes was developed for these sensors. The sensor’s design supports hot swapping and does not require restarting the system.
The method of contour mapping and source seeking in the radiation field by the UAS swarm equipped with gamma-ray detectors was developed. The method is used for low altitude applications where fixed wing UAS platforms are not suitable. The source-seeking and contour mapping algorithms were validated using a realistic radiation field with scattering and attenuation of gamma flux computed with MCNP code. The algorithm was implemented to map the radiation contours for multiple radiation sources and also a moving source. Moreover, it showed an effective way of cutting down the flight trajectories of three flying platforms by adaptively updating a swarm’s spin rate based on averaging the measured radiation data. The UAS swarm enables surveying an unknown, contaminated environment and mapping the area while locating the radiation sources, thus helping first responders to enhance situational awareness and to manage operations and safeguard personnel.
This work was supported by the Department of Energy Minority Serving Institution Partnership (MSIPP) managed by the Savannah River National Laboratory under SRNS contract MSIT00016.
Conflict of interest
The authors declare no conflicts of interest.