Spatial Deconvolution of Aerial Radiometric Survey and its application to the Fallout from a Radiological Dispersal Device
aa r X i v : . [ phy s i c s . d a t a - a n ] J a n Spatial Deconvolution of Aerial Radiometric Surveyand its application to the Falloutfrom a Radiological Dispersal Device
Laurel E. Sinclair, Richard Fortin
Canadian Hazards Information Service, Natural Resources Canada, Ottawa, Ontario,Canada
Abstract
Mapping radioactive contamination using aerial survey measurements is anarea under active investigation today. The radiometric aerial survey tech-nique has been extensively applied following reactor accidents and also wouldprovide a key tool for response to a malicious radiological or nuclear incident.Methods exist to calibrate the aerial survey system for quantification of theconcentration of natural radionuclides, which can provide guidance. However,these methods have anticipated a spatial distribution of the source which islarge in comparison to the survey altitude. In rapid emergency-responseaerial surveys of areas of safety concern, deposits of relatively small spatialextent may be expected. The activity of such spatially restricted hot spotsis underestimated using the traditional methods. We present here a spatialdeconvolution method which can recover some of the variation smoothed outby the averaging due to survey at altitude. We show that the method canrecover the true spatial distribution of concentration of a synthetic source.We then apply the method to real aerial survey data collected following deto-nation of a radiological dispersal device. The findings and implications of thedeconvolution are then discussed by reference to a groundbased truckbornesurvey over the same contamination.
Keywords: aerial, airborne, mobile survey, unfolding, deconvolution,inversion, MINUIT, MINOS
Full article appears in J. Environ. Radioact. 197 (2019) 39-47 . Introduction
Aerial radiometric survey is a mature field. Successful prospecting for eco-nomically viable ore deposits using the radiation signal from naturally occur-ing rocks stretches back decades [1]. Survey systems composed of large vol-umes of NaI(Tl) scintillator gamma-ray detectors, as much as 20 L, coupledwith georeferenced position sensors (now making use of the global naviga-tion satellite system (GNSS)), record gamma energy spectra versus position.This information is later processed to produce maps of natural potassium,uranium and thorium concentrations. Standards exist to guide practitionersin this area [2, 3] and vast regions of the earth have been covered [4]. Practi-tioners have also developed methods to correct for terrain variation in aerialsurvey [5, 6].The emphasis in aerial radiometric survey methods until recently hasbeen on development of techniques suitable for geologic sources, for whichthe simplification of the source as an infinite and uniform sheet is reasonablein comparison with the distance scales of the survey parameters (altitude,line spacing). The higher an aircraft flies, the more that far-away locationscontribute to the detected signal, relative to locations directly underneath theaircraft [7]. This can have the advantage of allowing for complete coveragein a more economical survey with wider line spacing. However, detectionsystems at higher altitude see a signal which is effectively averaged over alarger area of the ground. Anthropogenic signals such as those resulting incase of a reactor accident or malicious radiological dispersal could result inhot spots the concentration of which would be underestimated if averagedover a larger area.In this paper, we present a method to deconvolve an aerial radiometricsurvey for spatial smearing. This kind of problem, requiring inversion of aspatial distribution, is encountered frequently in geophysical surveying. Geo-physical spatial inversion problems are typically underdetermined, and oneway of dealing with this has been to select only those solutions which areclose to some preconceived model [8]. An approach to spatial deconvolutionof airborne spectrometric data which relies on an analytical model for theresponse function, and allows underdetermined problems, has been publishedpreviously [9]. A related method for spatial inversion to the approach pre-sented here, but using an iterative inversion and neglecting uncertainties, waspublished recently [10]. Other groups are taking a similar approach to thatadvocated here, but applied to the inversion of spectra rather than spatial2aps [11, 12].The method which will be presented here was applied to data obtained us-ing detectors aboard a manned helicopter. Nevertheless, it is prudent to men-tion the proliferation of work ongoing currently in aerial radiation detectionfrom unmanned aerial vehicle (UAV) systems. The use of a UAV platformfor aerial survey has facilitated the advance of aerial survey methodology. Agood review of recent publications can be found here [13].The method which will be presented here involves simply a) determiningthe influence of each independent region of the earth’s surface on the mea-sured spatial distribution and then b) optimizing the weight coefficients ofeach region of the surface to obtain the best reproduction of the measuredmap. The number of pixels of the solution can be chosen such that the prob-lem is not underdetermined. No potentially biasing prior assumption aboutthe underlying distribution is necessary. The method can handle complicateddetector geometries as the response functions are determined in simulation.The method could easily be extended to allow it to be applied when thereis significant terrain variation in the source such as would be the case in anurban environment. The method is stable under different starting conditions,and naturally allows for propagation of uncertainties from the measurementto the unfolded result.We demonstrate the application of the unfolding method by applying itfirst to a synthetic data set. This is compared with the known underlyingdistribution. We proceed to apply the unfolding method to a real aerialsurvey measurement acquired following detonation of a radiological dispersaldevice [14].This spatial deconvolution method has been presented previously at con-ferences by this group [15, 16], however, this is the first full write-up.
2. Methods
A measurement of surface activity concentration under the uniform infi-nite plane assumption may be denoted g MEAS ( x, y ) where x and y representeasting and northing in geographic Cartesian coordinates. We seek to deter-mine the true underlying surface radioactivity concentration, f ( x, y ). f ( x, y )is related to g MEAS ( x, y ) through g MEAS ( x, y ) = S [ f ( x, y )] , (1)3here S represents the effect of the measurement system on f ( x, y ).We divide space into N PAR pixels i , and using Monte Carlo simulation,generate uniform radioactivity distributions in each pixel, f i ( x, y ).The measurement system S consists of the detection system as well as theair and all other absorbing and scattering materials between the source andthe products of scattering, and the detectors. It is represented in the MonteCarlo simulation, and the emissions from the radioactive sources f i ( x, y ) aretransported through the system S to obtain the template responses g i ( x, y ),where g i ( x, y ) = S [ f i ( x, y )] . (2)We let g ( x, y ) = N PAR X i =1 w i g i ( x, y ) (3)and fit g ( x, y ) to g MEAS ( x, y ) using a χ minimization [17] to extract theweighting coefficients w i .To examine this χ function, let g MEAS j ( x, y ) represent the j th measure-ment of the activity concentration g MEAS ( x, y ). Then χ = N MEAS X j =1 g MEAS j ( x, y ) − P N PAR i =1 w i g i ( x, y ) e j (4)where there are N MEAS measurements g MEAS j ( x, y ) in the problem each withuncertainty e j .The estimator of f ( x, y ) is then the reconstructed radioactivity concen-tration distribution f REC ( x, y ), where f REC ( x, y ) = N PAR X i =1 w i f i ( x, y ) . (5)We choose to require the problem to be oversampled. That is, there iseverywhere a greater spatial density of measurements than of fit parametersand N MEAS > N
PAR . Then, provided that the uncertainties e j in the denom-inator of the χ function encompass all of the uncertainties of the problem,the minimum χ value will be approximately equal to the number of degreesof freedom of the problem. And the MINOS algorithm [17] can be used topropagate the N MEAS measurement uncertainties e j through the fit proce-dure to calculate the N PAR uncertainties δw i on the weighting parameters4 i . In practice, there are irreconcilable nonstochastic uncertainties affectingthe problem which must be included in the e j by application of a constantscaling factor to bring χ per degree of freedom to one before the fit uncer-tainties can be utilized. These are due to the statistical uncertainties in thetemplate responses g i ( x, y ), and the finite pixellization of the problem.Uncertainties for spatial deconvolution of fallout surveys can be expectedto be asymmetric owing to the boundary condition that the amount of de-position can not physically be less than zero. MINOS works by setting thepositive and negative uncertainty for each parameter to the amount the pa-rameter has to vary in each direction such that χ increases by one. Thus MI-NOS naturally allows for assymetric uncertainties and is particularly suitedto uncertainty analysis in the measurement of amount of radioactivity. The experimental methods to obtain the data to which we will apply theunfolding method have been described previously [14]. We will repeat onlythe most essential points here. Three RDD detonations were conducted dur-ing the trial. In the first, ∼
31 GBq of La-140 was dispersed explosively, withradioactive debris subsequently carried by wind as far as ∼ NaI(Tl) crystals mounted exteriorto a helicopter in a skid-mounted cargo expansion basket. GNSS antennaeand inertial navigation and altimetry systems were also installed in the bas-ket, to determine location. The system recorded a linearized 1024-channelgamma-ray energy spectrum over the domain 0 MeV to 3 MeV, tagged withthe location information, once per second. Post-acquisition, counts were se-lected from the spectra in an energy window of approximately four sigmain width around the 1.6 MeV La-140 photopeak. These count rates werecorrected for lag and dead time. Count rates were also corrected for varia-tions in flight altitude to the nominal flying height making use of the ShuttleRadar Topography Mission [18] digital terrain model adjusted to the localelevation using a GNSS base station. Backgrounds due to the radioactivityof the earth, the atmosphere, the aircraft and its occupants, and cosmic rays,were all subtracted. The count rates were all corrected for radioactive decay,back to the time of the blast. A coefficient to convert the measurementsfrom counts per second to kBq/m , assuming an infinite and uniform source,was obtained from experimentally validated Monte Carlo simulation. Finally,5easurements of radioactivity concentration in kBq/m for four aerial sur-veys, two conducted after the first blast and one conducted after each of thesubsequent two blasts, were presented.In this paper, we will discuss only the data recorded during the first aerialsurvey after the first blast. This survey was flown at a nominal 40 m flyingheight, with speed of 25 m/s and flight-line spacing of 50 m. Truckborne surveys were driven by criss-crossing the deposited falloutin an extemporaneous pattern following the first and third RDD blasts, re-stricting to the part of the fallout outside of a 500 m x 500 m fenced hotzone [19, 20]. The detection system was mounted in the bed of a pickuptruck and consisted of four 10.2 x 10.2 x 40.6 cm NaI(Tl) crystals orientedvertically in a self-shielding arrangement for azimuthal direction measure-ment. Truckborne data following the first RDD blast will be presented herefor comparison with the aerial survey data. The truckborne data has notundergone sufficient analysis for a full quantitative evaluation, but the shapewill nevertheless provide an interesting comparison for interpretation of theaerial survey data.
The sensitivity of the truckborne system to a La-140 contamination onthe ground was determined using experimentally validated Monte Carlo sim-ulation. In the simulations, the detector was placed with the centre of itssensitive volume at a height of 1.2 m above ground. The NaI(Tl) crystalsand their housing were represented in the simulation in their vertical arrange-ment. The ∼ R of a disc-shaped source centered be-neath the truckborne detector on the ground as determined by the simulation.The expression for the flux, Φ( R, H ), due to a surface activity concentration [m] R ) ] - m ⋅ / ( k B q - S en s i t i v i t y [ s EGSnrc
C(R,H) fit function asymptote
C(R,H) a) [m] R R e l a t i v e S en s i t i v i t y truckborne systemairborne system b) Figure 1: a) Sensitivity of the truckborne survey system to a disc-shaped distribution ofisotropic emitters of the La-140 gamma spectrum, as a function of disk radius. Black dotsshow EGSnrc prediction. Solid curve shows fit of the expression C ( R, H ) to the syntheticdata. The dashed line shows the asymptote of the fit curve and represents the sensitivityto a uniform and infinite sheet source. b) Comparison of the shapes of the sensitivitycurves for the aerial and truckborne survey systems. The dashed line shows the sensitivityto a disc source relative to the sensitivity to an infinite sheet as a function of disc radiusfor an aerial survey system at an altitude of 40 m [14]. The solid line shows the equivalentrelative sensitivity curve for the truckborne survey system. S at a point an elevation H above a disc-shaped source of radius R can bereadily calculated [7],Φ( R, H ) = S E ( µ a H ) − E ( µ a √ H + R )) , (6)where E is the exponential integral and µ a is the linear attenuation coefficientfor gamma rays in air. To determine the asymptotic sensitivity, we formed afunction for the detected counts as a function of the source radius, C ( R, H ) = ǫ Φ( R, H ) , (7)7nd fit the expression for C ( R, H ) to the synthetic data to obtain the detec-tion efficiency ǫ . The fit result is shown as the solid curve in Fig. 1a) and theasymptotic sensitivity, shown by the dashed line, is ∼
71 s − /(kBq/m ). Asmentioned, the uncertainty on this sensitivity is large due to the lack of de-tailed representation of the shielding material of the truck in the simulation.Nevertheless, the shape of the sensitivity curve is of value, as is the shapeof the profile of counts measured with the truckborne system as it traversedthe deposited plume. Fig. 1b) shows a comparison of the shapes of the sensitivity curves of theground-based and aerial systems. Despite the tall narrow shape of its detec-tors, which would tend to increase sensitivity to incoming radiation from thesides, for the truckborne system at “altitude” H = 1 . H = 40 m. This leadsto superior spatial precision in the results of truckborne survey. The radiation transport model EGSnrc [21, 22] was used to generatethe individual uniform pixel sources f i ( x, y ) and to propagate the generatedgamma rays and their interaction products through air and into the detectionvolume to create the template responses g i ( x, y ). For the solutions presentedherein, the simulation geometry represented the experimental setup duringthe first RDD trial [14]. The actual aerial survey system was described andshown in photographs in the previous publication and briefly reiterated inSect. 2.2. The model of that system in EGSnrc is shown in Fig. 2 a). Thesimulated gamma detection system included the two NaI(Tl) crystals, aswell as their aluminum cladding, and felt, foam and carbon fibre enclosure.The exterior-mounted basket containing the detectors was modelled in asimplified manner with 51 3 mm x 1.5 mm bars of aluminum, representingthe basket strands, running the length of the basket in a semicircle aroundits long axis. Dead materials due to the photo-multiplier tube readout ofthe crystals and associated electronics, as well as the altimeter and GNSSreceivers, were modelled as simple blocks of metal of the appropriate overallmass. The ground was represented as perfectly flat with the detection systemat a height of 40 m. The model was validated experimentally using data8rom point sources of known emission rate placed at known locations aroundthe actual detector. The uncertainty associated with approximation in therepresentation of the measurement system in the model was evaluated byvariation of the arrangement of dead material in the model, by variation of thedetector’s position, altitude and attitude and by variation of the detector’senergy resolution within reasonable limits. This uncertainty was evaluatedto be about 12% on the activity concentration and it was included the overallsystematic uncertainty quoted in the publication [14].a) Easting [m] N o r t h i ng [ m ] % -6 × b) Figure 2: a) The aerial survey system as modelled in EGSnrc. The two10.2 x 10.2 x 40.6 cm NaI(Tl) crystals are represented in their housing, shown in purple,with steel blocks at the ends representing the dead material of the PMTs and readoutelectronics. The two crystals are mounted lengthwise with their centres 152.4 cm aparton an aluminium plate which is represented with mass and dimensions according to theengineering drawing. Auxilliary instrumentation is represented by an aluminium block inthe centre of the basket, of the summed instrument mass. This dead material is shown ingreen in the figure. The aluminum plate on which the detectors and auxiliary equipmentis mounted is itself attached to a basket which is mounted to the skids on the outsideof the helicopter. Dead material of the basket is represented by 51 thin aluminum barsrunning the length of the basket, in a semicircle around the basket long axis, shown ingold. Uncertainty on the sensitivities due to misrepresentation of the system in the modelhas been estimated to be approximately 12%. b) Response function, g ( x, y ) to the pixelsource f ( x, y ), normalized to the number of generated events, where the x -axis showsEasting and the y -axis shows Northing. The true spatial extent of pixel 53 is indicated bythe black square. The simulated pixel sources, f i ( x, y ), consisted of uniform distributions ofisotropic emitters of the La-140 spectrum of gamma rays, including all emis-9ion energies above 0.05% emission probability [23]. Gammas were emittedinto the full 4 π solid angle. Each pixel source f i ( x, y ) was square, and 50 mon a side. Note that with the survey parameters mentioned previously, wehave one spectrum recorded approximately every 25 m x 50 m in the realdata. Thus, the number of measurements in the data is higher than thenumber of fit parameters in the simulation, the problem is over-determined,and we can expect reasonably rapid convergence of the method to a solutionwhich is stable under different starting conditions.The entire region to be unfolded measured 1.5 km x 1.5 km. To speed upconvergence of the fitting, we chose to parallelize the computation, breakingthe problem up into individual tiles, each 500 m x 500 m. Only eight of thesetiles were necessary to cover the area over which the radioactivity actuallysettled. To knit the tiles together, we extended the fit into a border of pixelsources 50 m wide, surrounding each tile. Thus the fit area corresponding toeach tile was actually 600 m x 600 m. If, after background subtractions, fewerthan one count corresponding to a La-140 energy deposit was measured inthe detection system in the region of space corresponding to one fit parameter w i , then that fit parameter was assigned a value of zero and left out of theminimization. Thus, 144 or fewer fit parameters w i were determined fromthe inversion of each tile. To merge the tiles after the individual inversions,the weighting parameters of the border pixels were simply ignored, and thecentral 500 m x 500 m areas of the tiles were placed next to each other.The simulation included one detection system for each of the 144 pixelsources, centered laterally within the pixel, at 40 m height. Here we have usedmultiple detection systems at different places at one time to represent the realworld in which one detection system was in different places at different times.Given the 40 m height of the detection systems above the source, the 50 mspacing between them and the requirement that the deposited energy lie inthe highest-energy photopeak of the source, the error in this approximation,which would come from a single emitted gamma ray depositing energy in twodifferent detection systems, is negligible.The volume of the air in which the gamma rays and electrons were trackedin each tile extended 1.2 km in easting (or x ) and northing (or y ), and 600 mup. A 5 m-thick concrete slab underneath the emitters and filling the lateraldimensions of the simulated volume, represented the earth.The simulation included all physical interactions of the emitted gammasand of the gammas, electrons and positrons resulting from those interactions.Scattering and absorption in the air and “earth” of the simulated volume,10nd in the dead material of the basket system and housing of the NaI(Tl)crystals, was included. All processes leading to energy deposit within theNaI(Tl) crystal from either the direct gamma rays, or the products of scat-tering, were included. Energy deposits in the NaI(Tl) were then smeared tocreate simulated spectra, according to the energy resolutions of the crystalsdetermined during the experiment.Fig. 2 b) shows the response, g ( x, y ), of the 144 detection systems ofone tile, as a percentage of the number of events generated in a single oneof the pixel sources, f ( x, y ), where this happens to correspond to the pixelnumbered “53”. The extent of the source is indicated by the black square.Note how the measured response extends much more broadly in space. Thisis the spatial smearing which will be corrected by the unfolding.
3. Results
We begin by applying the spatial inversion method to a synthetic dataset for which we know the underlying distribution f ( x, y ). An annular regionconsisting of the area between two circles of radius 100 m and 250 m, centeredat (550 m, 550 m) was uniformly populated with 10 · emitters of theLa-140 gamma sectrum. Considering the branching ratios for the La-140gamma emissions, this amounts to a source concentration of 28.4 kBq/m .The annular region is indicated by the black outlines in Fig.’s 3 a) and b).Generation of the synthetic dataset makes use of the identical detectorsimulations as used in generating the template responses g i ( x, y ) as describedin Section 2.4, however the detection systems were more narrowly spacedin the synthetic dataset. Detection systems for the synthetic dataset werelocated every 20 m in x and y such that there were 225 detection systems intotal in each 600 m x 600 m tile.The template sources f i ( x, y ) and the template responses g i ( x, y ) utilizedin the inversion are the same as will be used for the real data and are asdescribed in Section 2.4. Thus, the synthetic data measurement density of20 m x 20 m exceeds the density of the parametrization, 25 m x 25 m, andthe problem is overdetermined, as desired.Fig. 3 a) shows the synthetic aerial survey measurement. The result isbroader than the underlying true source distribution. Contamination appearsto extend outside of the known true underlying borders. The central area ofthe annulus appears to be filled with significant contamination. The average11 asting [m]
200 400 600 800 1000 1200 1400 N o r t h i ng [ m ] ] S u r f a c e a c t i v i t y c on c en t r a t i on [ k B q / m a) Easting [m] N o r t h i ng [ m ] ] S u r f a c e a c t i v i t y c on c en t r a t i on [ k B q / m b) Figure 3: Synthetic aerial survey and fit results, where the x axis is Easting, and the y axis is Northing. The source was of concentration 28.4 kBq/m and annular with innerradius 100 m and outer radius 250 m, centered at (550 m, 550 m), as indicated by thearea between the black circles. a) Synthetic aerial survey measurement result. b) Resultof fit of template responses g i ( x, y ) to the synthetic aerial survey. concentration of contamination reconstructed within the annular region islower than the known true concentration.Fig. 3 b) shows the result of the fit of the template measured activitydistributions to the synthetic aerial survey measurement. The colour scaleused in Fig. 3 b) is the same as that of Fig. 3 a). The first observation to noteis that the tile knitting procedure apparently works well. The synthetic datastretches over six of the overlapping 600 m x 600 m simulation tiles. Afterknitting of the inverted result into adjacent non-overlapping 500 m x 500 mtiles, there is seamless matching of the reconstructed activity concentrationat the tile borders. Also to note is that within the limitations of the somewhatcoarse pixellization of the problem, the survey is well reproduced by the fit.The tendency of the measurement to extend outside of the bounds of the truesource is reproduced, as is the tendency of the measurement to underestimatethe magnitude of the concentration within the source boundary.The reconstruction of the true underlying source distribution for the syn-thetic data is presented in Fig. 4. As shown in Fig. 4 a), the inversionprocess results in a reconstructed source distribution which is higher in mag-nitude and closer to the true activity concentration than the initial surveymeasurement. The boundaries of the actual source are better reproducedafter inversion, in particular the absence of radioactivity in the centre of the12
200 400 600 800 1000 1200 1400
Easting [m] N o r t h i ng [ m ] ] S u r f a c e a c t i v i t y c on c en t r a t i on [ k B q / m a) Easting [m] N o r t h i ng [ m ] ] S u r f a c e a c t i v i t y c on c en t r a t i on [ k B q / m b) Easting [m] N o r t h i ng [ m ] ] S u r f a c e a c t i v i t y c on c en t r a t i on [ k B q / m c) Figure 4: Spatially deconvolved synthetic aerial survey result, x axis is Easting and y axis is Northing. Black circles indicate the true annular source distribution. a) Thespatially deconvolved measurement. b) Positive statistical uncertainty on the spatiallydeconvolved measurement. c) Negative statistical uncertainty on the spatially deconvolvedmeasurement. annulus is recovered.A major advantage of the spatial deconvolution method presented here isthat statistical uncertainties affecting the measurement may be propagatedthrough the minimization procedure by the MINOS algorithm as describedin Sect. 2.1. A map showing the one-sigma positive statistical uncertainty onthe reconstructed surface activity concentration is shown in Fig. 4 b) and thecorresponding negative statistical uncertainty is shown in Fig. 4 c) where thesame colour scale is used for the uncertainties as for the measurement shownin Fig. 4 a). Considering the uncertainties, the ability of the method to recon-struct the true underlying activity concentration is good. The reconstructedactivity concentration magnitude is in agreement with the known true activ-13ty concentration within uncertainties in most places. For example, near theinner boundary of the annulus where the reconstructed concentration is lowcompared to the known true concentration, further negative movement of theresult is not allowed by the negative uncertainty. The positive uncertaintyexceeds the negative uncertainty in magnitude, and does allow for a positivemovement of the reconstructed value. The aerial survey measurement from the first RDD trial is shown inFig. 5 a). This result has been published previously [14] and the meth- ] S u r f a c e a c t i v i t y c on c en t r a t i on [ k B q / m Easting [m] N o r t h i ng [ m ] ] S u r f a c e a c t i v i t y c on c en t r a t i on [ k B q / m a) ] S u r f a c e a c t i v i t y c on c en t r a t i on [ k B q / m Easting [m] N o r t h i ng [ m ] b) Figure 5: a) Aerial survey measurement of the distribution of fallout following detonationof the radiological dispersal device. b) Result of fit of template histograms to the aerialsurvey measurement. ods to arrive at the result were recapitulated here in Sect. 2.2. We observean area of activity concentration exceeding 100 kBq/m near ground zero ofthe detonation, with a long deposited plume of maximum width of 300 mto 400 m, and significant measured radioactivity extending over a distanceof about 2 km. The total systematic uncertainty affecting this measurementwas determined to be around 12% and the statistical uncertainty peaks at6 kBq/m .Fig. 5 b) shows the result of the χ fit of the weighting coefficients of thetemplate response functions to the measurement of Fig. 5 a). The colourscales of Fig.’s 5 a) and b) have been chosen to be equal. (This colour scalehas in fact been set to the optimal colour scale for the data from a truckborne14urvey which will be presented in the upcoming Sect. 3.2.2.) The featuresof the measurement are broadly well-reproduced by the fit, considering thepixellization of the reconstruction. In particular, the magnitude, width andextent of the contamination are well represented in the result of the fit.The underlying distribution of pixel sources which gives rise to the fitresult of Fig. 5 b) is presented in Fig. 6 a). This is the reconstructed distri- ] S u r f a c e a c t i v i t y c on c en t r a t i on [ k B q / m Easting [m] N o r t h i ng [ m ] a) ] S u r f a c e a c t i v i t y c on c en t r a t i on [ k B q / m Easting [m] N o r t h i ng [ m ] b) ] S u r f a c e a c t i v i t y c on c en t r a t i on [ k B q / m Easting [m] N o r t h i ng [ m ] c) Figure 6: a) Spatially deconvolved aerial survey measurement of fallout following radio-logical dispersal device blast. b) Positive statistical uncertainty on spatially deconvolvedmeasurement. c) Negative statistical uncertainty on spatially deconvolved measurement. bution of the surface activity concentration following spatial deconvolution.We find that the width of the deposited plume is now much smaller thanthe original undeconvolved measurement, around 50 m. Correspondingly,the peak concentration is higher, over 400 kBq/m . Note that the widthof the deposited plume is small with respect to the altitude and line spac-ing of the survey. This accounts for its significant overestimation when the15infinite and uniform” approximation was used to obtain a concentrationmeasurement from the measured counts as shown in Fig. 5 a) and [14]. Thelength of the deposition is however much larger than the survey parameters,so in this dimension the “infinite and uniform” sheet approximation is not sobad and the original length measurement is not much altered by the spatialdeconvolution.The positive and negative statistical uncertainties on the spatially decon-volved deposited fallout map are shown in Fig.’s 6 b) and c) respectively.Note that the statistical uncertainties affecting the measurement are verysmall, and on the colour scale of Fig. 6 a) would be difficult to see, so thecolour scale in Fig.’s 6 b) and c) is chosen to optimize the representationof the information in Fig. 6 c). The uncertainties reveal important featuresof the measurement and its inversion. The positive uncertainty indicates aregion extending away from the measured deposited plume axis in which apositive quantity of activity is permitted, however at a very low amount ofbetween 5 kBq/m and 10 kBq/m . The negative uncertainty shows that themeasurement of the presence and distribution of radioactivity is significantlyabove zero.The MINOS error propagation includes only the stochastic uncertaintieson the measurement. There are additional uncertainties which are systematicand arise from approximations in the representation of the system in thesimulation. These include mis-representation of the position, particularly thealtitude; the attitude (yaw, pitch and roll); the amount of shielding materialin the basket containing the detectors; the energy resolution and the airdensity. The systematic uncertainty on the (undeconvolved) radioactivityconcentration distribution was determined to be about 12% by variation ofthese parameters within reasonable limits [14].For the spatially deconvolved measurement, some of the systematic un-certainties must be re-examined as they can be expected to have an effecton the shape of the reconstructed fallout distribution, as well as its overallmagnitude. These are the systematic uncertainties associated with the mea-surement of altitude and pitch angle. It is also interesting to examine the ef-fect of the measurement of yaw angle on the spatially inverted measurementas it can have no effect at all on the original undeconvolved measurementwhich used the infinite sheet approximation for the source and therefore thedetector response was invariant under changes of yaw.The inertial navigation system determined the yaw angle during the mea-surement to be around -30 ◦ . The spatial inversion was conducted twice. For16he central value of the inversion as presented in Fig. 6 a) the helicoptersystems in the template histograms were assigned a yaw of -30 ◦ to matchthe data. To allow for changes in yaw during flight, the regression was re-peated with yaw set maximally different at 60 ◦ . Pitch was varied from thenominal 0 ◦ to -10 ◦ according to the maximum deviation of pitch recordedby the inertial navigation system while on line during one sortie. Altitudewas varied from nominal by 1 m to account for approximately one sigmaof variability in height determination. These variations did not significantlyalter the measurement of length and width of the deposited fallout. Addedin quadrature, and considering that some of the variation was already in-cluded in the original systematic uncertainty associated with the sensitivity,the deviations do not yield a significant additional systematic uncertainty.Although not a significant additional source of uncertainty for the measure-ments presented here, these sources of uncertainty are worth discussing forthe benefit of researchers following this approach under different operatingconditions. In Fig. 7a) the truckborne survey result which followed the first blast isshown overlaid on the undeconvolved aerial survey result from the same blaston the same colour scale. (The colour scale is optimized to show the variationin the truckborne survey result.) The truckborne survey reports much highersurface activity oncentrations than the aerial survey and the fallout appearsto be narrower in width.Fig. 7b) shows the same truckborne survey result this time overlaid onthe spatially deconvolved aerial survey measurement. (The colour scale isthe same as that used in Fig. 7a).) Here, both the width and magnitude ofthe concentration are in better agreement.The aerial survey maps were sampled at the locations of the truck andthese sampled activity concentration values are presented in Fig. 7c) (top) asa function of the distance from the start of the truck-driven sortie. Again, itis clear that the spatially deconvolved result is generally higher at maximummagnitude and more narrow than the undeconvolved aerial survey result.The bottom part of Fig. 7c) shows the root-mean-square deviation (RMSD)of the curves calculated by sampling each profile at regular intervals from themaximum concentration down to the point at which the concentration fallsbelow 10% of the maximum. This plot shows that the RMSD width of the17eposition according to the undeconvolved survey is typically around 90 mwhile the width of the deposition according to the truckborne survey and thespatially deconvolved aerial survey tends to be significantly narrower, closerto 30 m. Spatially deconvolved aerial survey thus recovers the narrowness ofthe fallout to approximately the same spectral precision as the truckbornesurvey.Truckborne data is shown for the purpose of shape comparison only anddoes not include detailed error analysis. In any event, the truckborne systemhas its own finite area of sensitivity largely caused by its “altitude” of justover one metre, causing smearing of the measured spatial distribution. Acontact measurement of the deposited radioactivity can be expected to beeven more concentrated in places [24].
4. Discussion and Conclusions
Radiometric survey would be performed to map fallout following a reactoraccident or following a malicious release. To cover a large area quickly, thesurveys are initially performed using manned aircraft at some significantaltitude H . If there is spatial variability in the fallout at distance scalesmuch less than H , then the map result of the survey can underestimate thequantitity of radioactivity on the ground in places.We have presented here a method to deconvolve an aerial survey map forthe spatial smearing caused by measurement at altitude, at least to the ex-tent permitted by the sampling density as determined by the aircraft speedand line spacing. Performed on synthetic data, the deconvolution methodreturns a distribution which is consistent with the true underlying distribu-tion within measurement uncertainty. The deconvolved distribution is morenarrowly distributed, and shows regions of locally higher radioactivity con-centration than the initial undeconvolved measurement. Performed on realdata acquired following detonation of a radiological dispersal device, themethod produces a distribution which is narrower and shows radioactivityconcentration as much as four times that of the original measurement.The method can unfold a distribution for smearing effects up to a reso-lution permittable by the sampling frequency of the original measurement.The method allows for propagation of stochastic measurement uncertaintiesthrough the unfolding to obtain the measurement uncertainties on the fitparameters. The method relies on application of the MINUIT and MINOSalgorithms well known in the field of particle physics. What is perhaps not18ell known is that these algorithms can tolerate operating with hundreds ofindependent fit parameters, converging to a stable solution in a reasonableamount of time from an arbitrary starting distribution. Our current needwas to develop a method to extract the greatest possible information froma set of aerial surveys performed to improve scientific understanding of thebehaviour of radiological dispersal devices. The method, however, is alsoapplicable to unfolding of any smeared distribution of any dimensionality. Itcould find application in other fields.The result of the unfolding is limited in spatial resolution by the require-ment that the density of pixellization of the answer not exceed the densityof measurements as determined by the aircraft speed and line spacing. Nev-ertheless, aerial survey practitioners should be aware that there is improvedinformation about the spatial distribution of the radioactivity contained intheir aerial survey map that can be extracted provided good knowledge ofthe response function of the system is available. The achieved spatial reso-lution for the particular aerial survey following RDD detonation presentedhere approximately matched that obtained during truckborne survey over thesame deposition (while providing complete coverage). Contact measurements( H = 0 m) can be expected to reveal even greater local spatial variationsthan the truckborne data. Still, the truckborne survey “height” of about1.2 m provides a salient benchmark for spatial resolution as this is close tothe average height of an adult human. Should humans be required to entera possibly contaminated area guided by the results of aerial survey alone,a spatially deconvolved aerial survey map could provide a better predictorof the activity concentrations they will encounter than the undeconvolvedmeasurement. Acknowledgements
The authors gratefully acknowledge the leadership of the RDD field trialsunder L. Erhardt, and helpful comments on the analysis from H.C.J. Seyw-erd, P.R.B. Saull and F.A, Marshall. Funding for this project was providedthrough Canada’s Chemical, Biological, Radiological-Nuclear and ExplosivesResearch and Technology Initiative. This report is NRCan Contribution20180112. 19 eferences [1] R. L. Grasty, Uranium measurement by airborne gamma-ray spectrom-etry, Geophysics 40 (1975) 503.[2] Airborne gamma ray spectrometer surveying, IAEA Technical ReportsSeries 323 (1991).[3] R. L. Grasty, Environmental monitoring by airborne gamma ray spec-trometry, experience at the Geological Survey of Canada, Applicationof uranium exploration data and techniques in environmental studies,IAEA-TECDOC-827 (1995).[4] Radioelement Mapping, IAEA Nuclear Energy Series No. NF-T-1.3(2010).[5] G. Schwarz, E. Klingele, L. Rybach, How to handle rugged topographyin airborne gamma-ray spectrometry surveys, First Break 10 (1992) 11– 17.[6] A. Ishizaki, Y. Sanada, M. Ishida, M. Munakata, Application of topo-graphical source model for air dose rates conversions in aerial radiationmonitoring, J. Environ. Radioactiv. 180 (2017) 82 – 89.[7] L. V. King, Absorption problems in radioactivity, Phil. Mag. 23 (1912)242–250.[8] R. L. ”Parker, ”Geophysical Inverse Theory”, ”Princeton UniversityPress”, ”41 William St., Princeton, New Jersey 08540, USA”, ”2015”.[9] B. Minty, R. Brodie, The 3D inversion of airborne gamma-ray spectro-metric data, Exploration Geophysics 47 (2016) 150 – 157.[10] R. D. Penny, T. M. Crowley, B. M. Gardner, M. J. Mandell, Y. Guo,E. B. Haas, D. J. Knize, R. A. Kuharski, D. Ranta, R. Shyffer, S. Labov,K. Nelson, B. Seilhan, J. D. Valentine, Improved radiological/nuclearsource localization in variable NORM background: An MLEM approachwith segmentation data, Nucl. Instr. and Meth. A 784 (2015) 319 – 325.[11] J. Mattingly, D. J. Mitchell, A framework for the solution of inverseradiation transport problems, IEEE Trans. Nucl. Sci. 57 (6) (2010) 3734– 3743. 2012] E. M. A. Hussein, The physical and mathematical aspects of inverseproblems in radiation detection and applications, Appl. Radiat. Isot. 70(2012) 1131 – 1135.[13] D. Connor, P. G. Martin, T. B. Scott, Airborne radiation mapping:overview and application of current and future aerial systems, Int. J.Remote Sens. 37 (2016) 5953 – 5987.[14] L. Sinclair, R. Fortin, J. Buckle, M. Coyle, R. Van Brabant, B. Harvey,H. Seywerd, M. McCurdy, Aerial Mobile Radiation Survey FollowingDetonation of a Radiological Dispersal Device, Health Phys. 110 (2016)458 – 470.[15] L. Sinclair, R. Fortin, Spatial Deconvolution of Aerial Radiometric Survey Results,CTBT: Science and Technology Conference, Presentation T3.2-P20,Vienna, Austria, accessed: March 23, 2018 (2015).URL [16] L. E. Sinclair, F. A. Marshall, R. Fortin, Reconstruction of the spatialdistribution of radioactive contamination from aerial survey and from astationary array of directional detectors, in: 2015 IEEE Nuclear ScienceSymposium and Medical Imaging Conference (NSS/MIC), 2015, pp. 1–4. doi:10.1109/NSSMIC.2015.7581989 .[17] F. James, M. Roos, Minuit: A system for function minimization andanalysis of the parameter errors and correlations, Comput. Phys. Com-mun. 10 (1975) 343–367. doi:10.1016/0010-4655(75)90039-9 .[18] T. G. Farr, P. A. Rosen, E. Caro, R. Crippen, R. Duren, S. Hens-ley, M. Kobrick, M. Paller, E. Rodriguez, L. Roth, D. Seal, S. Shaffer,J. Shimada, J. Umland, M. Werner, M. Oskin, D. Burbank, D. Alsdorf,The shuttle radar topography mission, Rev. Geophys. 45 (2007) –.URL http://dx.doi.org/10.1029/2005RG000183 [19] A. Green, L. Erhardt, L. Lebel, M. Duke, T. Jones, D. White, D. Quayle,Overview of the Full-scale Radiological Dispersal Device Field Trials,Health Phys. 110 (2016) 403 – 417.[20] F. Marshall, Reconstruction of the Spatial Distribution of Surface Ac-tivity Concentration for an In-Situ, Gamma-Ray, Truck-borne Survey ,Master’s thesis, Carleton University, Canada (2014).2121] I. Kawrakow, E. Mainegra-Hing, D. W. O. Rogers, F. Tessier, B. R. B.Walters, The EGSnrc Code System, NRC Report PIRS-701 (2011).[22] I. Kawrakow, E. Mainegra-Hing, F. Tessier, B. R. B. Walters, TheEGSnrc C++ class library, NRC Report PIRS-898 (2011).[23] M.-M. B´e, V. Chist´e, C. Dulieu, E. Browne, V. Chechev, N. Kuzmenko,R. Helmer, A. Nichols, E. Sch¨onfeld, R. Dersch, Table of Radionu-clides, Vol. 1 of Monographie BIPM-5, Bureau International des Poidset Mesures, Pavillon de Breteuil, F-92310 S`evres, France, 2004.[24] L. Erhardt, L. Lebel, E. Korpach, R. Berg, E. Inrig, I. Watson, C. Liu,C. Gilhuly, D. Quayle, Deposition Measurements from the Full-ScaleRadiological Dispersal Device Field Trials, Health Phys. 110 (2016) 442– 457. 22
200 400 600 800 1000 1200 1400