Modeling Aerial Gamma-Ray Backgrounds using Non-negative Matrix Factorization
M. S. Bandstra, T. H. Y. Joshi, K. J. Bilton, A. Zoglauer, B. J. Quiter
IIEEE TRANSACTIONS ON NUCLEAR SCIENCE 1
Modeling Aerial Gamma-Ray Backgrounds usingNon-negative Matrix Factorization
M. S. Bandstra, T. H. Y. Joshi, K. J. Bilton, A. Zoglauer, and B. J. Quiter
Abstract —Airborne gamma-ray surveys are useful for manyapplications, ranging from geology and mining to public healthand nuclear security. In all these contexts, the ability to de-compose a measured spectrum into a linear combination ofbackground source terms can provide useful insights into thedata and lead to improvements over techniques that use spectralenergy windows. Multiple methods for the linear decompositionof spectra exist but are subject to various drawbacks, such as al-lowing negative photon fluxes or requiring detailed Monte Carlomodeling. We propose using Non-negative Matrix Factorization(NMF) as a data-driven approach to spectral decomposition.Using aerial surveys that include flights over water, we demon-strate that the mathematical approach of NMF finds physicallyrelevant structure in aerial gamma-ray background, namely thatmeasured spectra can be expressed as the sum of nearby terres-trial emission, distant terrestrial emission, and radon and cosmicemission. These NMF background components are comparedto the background components obtained using Noise-AdjustedSingular Value Decomposition (NASVD), which contain negativephoton fluxes and thus do not represent emission spectra in asstraightforward a way. Finally, we comment on potential areasof research that are enabled by NMF decompositions, such asnew approaches to spectral anomaly detection and data fusion.
I. I
NTRODUCTION A ERIAL surveys are used to measure gamma-ray emis-sion over large areas, with applications in geologicalexploration, radiation protection, mapping distributed contam-ination, and finding radiological sources outside of regulatorycontrol [1], [2], [3], [4]. Commonly used analysis techniquesutilize the counts in certain spectral energy windows to mea-sure and map background levels or to detect sources, and theseapproaches generally trade statistics for specificity to particularradioactive isotopes [5], [6], [7]. Other techniques have beendeveloped to leverage the information contained across theentire spectrum, thus gaining statistics but potentially los-ing specificity [8], [9], [10], [11]. Though potentially morepowerful than windowed approaches, these techniques requiredetailed modeling of the different background components inorder to compensate for the loss of specificity. Here we willreview current full-spectrum approaches and present a newapproach using Non-negative Matrix Factorization (NMF) that,
M. S. Bandstra, T. H. Y. Joshi, and B. J. Quiter are with the NuclearScience Division at Lawrence Berkeley National Laboratory, Berkeley, CA94720 USA e-mail: [email protected]. J. Bilton is with the Department of Nuclear Engineering at the Universityof California, Berkeley, Berkeley, CA 94720 USA.A. Zoglauer is with the Berkeley Institute for Data Science and the SpaceSciences Laboratory at the University of California, Berkeley, Berkeley, CA94720 USA.Manuscript received —. Revised —. due to the reformulation of the problem, provides new insightsinto airborne survey measurements.
A. Aerial gamma-ray background
The natural backgrounds encountered by airborne systemsare a combination of terrestrial sources and non-terrestrialsources (e.g., [12]). The terrestrial sources are typically a com-bination of K, U/ U decay series isotopes, and
Thdecay series isotopes, collectively denoted KUT backgrounds.Emission from terrestrial sources can be classified into twotypes: direct emission, consisting of unscattered and scatteredphotons that are incident on the detector from below, and skyshine , consisting of photons that have scattered in the airabove the detector that are incident from above [10]. Thenon-terrestrial backgrounds are radon gas and its radioactiveprogeny, the effects of cosmic rays, and KUT backgroundsfrom the aircraft itself. The radon component consists pri-marily of
Rn progeny, and the cosmic component contains keV emission from cosmogenic positron production anda cosmogenic power-law continuum [13], [14].An illustrative study of the full-spectrum background com-ponents for a stationary ground-based detector was undertakenin refs. [13], [14]. In that study the authors empiricallyseparate terrestrial from non-terrestrial backgrounds, isolateskyshine from other terrestrial emission, and generate modelsthat correctly reproduce each component. Fully modeling thebackground required Monte Carlo simulations with terrestrialemission up to at least m away from the detector and with m of atmosphere overhead to produce the correct amountof skyshine, showing the complexity of the problem for evena static detector on the ground. Other ground-based studieshave shown the influence of skyshine at low energies even fordistant sources [15], [16].
B. Background models
Before going any further, we note that the term spectrummodel will carry a dual meaning in this work. In the firstsense, a spectrum model is a Monte Carlo physics simulationthat reproduces the components of a measured spectrum. Inthe second sense, a spectrum model is any mathematicaldecomposition of a measured spectrum, regardless of physicalmeaning of the components. Both models provide a numericaldescription of the data, with the former providing a physicalunderstanding and the latter providing a mathematical under-standing. An ideal method would include aspects of both kindsof models in one physical and mathematical description. c (cid:13) a r X i v : . [ phy s i c s . d a t a - a n ] F e b EEE TRANSACTIONS ON NUCLEAR SCIENCE 2
For a mathematical model of gamma-ray background to beconsistent with physics (e.g., the additive nature of photons),the measured spectra should be expressed as a linear sum ofspectral components. For this work, we begin with a measureddataset X , which is an n × d matrix where each of the n rows is a measured spectrum with d bins. We consider lineardecompositions of X of the form X ≈ ˆX = AV , (1)where each row of the n × d matrix ˆX is a model-estimatedspectrum, A is an n × k matrix of weights, V is a k × d matrix whose rows are the component spectra, and k ≤ d isthe number of components. In nearly all cases, equation (1) isan ill-posed problem with many potential solutions that dependon the choice of k and desired constraints on A and V . Herewe will review some of the approaches for finding solutions,along with their advantages and disadvantages. C. Full Spectrum Analysis
Full Spectrum Analysis (FSA) of gamma-ray spectra is oneapproach that has been studied for aerial and borehole detectorsystems [8], [9], [10], [17], [18]. FSA aims to express a mea-sured spectrum as a linear combination of component spectrathat are predetermined through simulations and calibrationmeasurements. For example, component spectra for the KUTbackgrounds have been developed to determine the isotopiccompositions of the ground beneath an aircraft [18]. The k spectral components make up the rows of V , and the weights A are typically calculated by minimizing the sum of the squaredifferences between X and ˆX , but since least squares canadmit non-physical (i.e., negative) solutions to equation (1),using non-negative least squares has been suggested to enforcethe non-negativity of A [19]. Another desirable propertyfor a spectrum model is to preserve the number of countsin each spectrum, which least squares minimization alonedoes not enforce, but this constraint can be included in theminimization [19].One challenge of FSA is that the decomposition can bepoor if the component spectra are not correct, so care mustbe taken to ensure their accuracy. The component spectra areusually determined by using a combination of Monte Carlosimulations and calibration measurements, for example bymeasuring slabs of material with known KUT concentrationsunderneath wooden boards to simulate direct terrestrial emis-sion (e.g., [20]). However, such measurements are unable toreproduce the effects of skyshine, which is a major terrestrialcontribution to airborne spectra below keV [10]. D. Noise-Adjusted Singular Value Decomposition
Other full-spectrum modeling techniques have focused onusing mathematical methods to discover and utilize structurefound within measured data without requiring a physicalunderstanding. One commonly used method is Noise-AdjustedSingular Value Decomposition (NASVD), which was devel-oped to remove statistical noise from aerial gamma-ray mea-surements by attempting to separate true spectral variationsfrom statistical noise [21], [22], [23]. NASVD improved upon similar work using Principal Component Analysis (PCA) [20].NASVD has been used to map low levels of contamina-tion [11] and has recently been applied to the problem of esti-mating real-time backgrounds for aerial surveys [24]. NASVDdecomposes measured spectra into principal components andthen keeps only those components which represent true spec-tral variability. Although NASVD is typically used to smoothspectra before further processing, isotope-specific signaturesin the first few components have led some to comment on thepossible associations between NASVD components and KUTand radon backgrounds (e.g., [22]). However, unlike FSA, allbut the first component contain negative values, so this methodpermits solutions to equation (1) that have negative photonfluxes, which are non-physical.
E. Non-negative Matrix Factorization
Recently Non-negative Matrix Factorization (NMF) [25],[26] was introduced as a technique for performing dimension-ality reduction of gamma-ray background data [27] and forgamma-ray source identification [28]. NMF has already beenused for spectral deconvolution in other fields of spectroscopy(e.g., [29], [30], [31]). NMF is a method of solving equa-tion (1) such that both A and V are constrained to benon-negative. NMF can be thought of as a bridge betweenFSA and NASVD — like NASVD, the NMF componentsare determined from the measurements themselves, and likeFSA, the NMF components are compatible with physics andthus are capable (although not guaranteed) of describing actualbackground emissions. Indeed recent analysis of data from avehicle-borne system has found correlations between NMF-derived background components and features derived frompanoramic images of the scene around the vehicle [32]. Otherwork has shown NMF to have advantages over other methodsfor spectral anomaly detection and identification [27].In this work NMF is applied to aerial survey data and theresulting components are interpreted as known backgroundsources, thus connecting a purely mathematical model ofthe spectra with a physical model. We start in Section IIby analyzing a survey containing passes over a land-waterinterface at low altitude to demonstrate the basic method.NMF components are compared to a full Monte Carlo modelin Section II-C and Section II-D, and then compared withNASVD components in Section II-E. More complex datasetsare analyzed in Section III, showing the potentially wide appli-cability of the method. Finally, the implications for mapping,anomaly detection, and data fusion using NMF are discussedin Section IV.II. A NALYSIS OF A LAND - WATER INTERFACE
The data analyzed here were taken on 14 February 2012 atLake Mohave on the Nevada-Arizona border. The detector sys-tem was flown on a helicopter and consisted of four RadiationSolutions Inc. RSX-3 units, totaling twelve × × -inchNaI(Tl) detectors [33]. Photon events were recorded by thedata acquisition system in keV bins from to keV,and for this analysis an energy cut of − keV wasused. Events were rebinned into d = 200 non-uniform bins EEE TRANSACTIONS ON NUCLEAR SCIENCE 3 -114.700 -114.690 -114.680 -114.670 -114.660Longitude+35.4350+35.4400+35.4450+35.4500+35.4550+35.4600+35.4650 L a t i t u d e Figure 1. The helicopter flight path of the Lake Mohave dataset showingrepeated land-water crossings. The inset shows a magnified view of theshoreline where the land-water crossings were made. (Map imagery: Google.) with bin widths approximately proportional to the squareroot of energy in order to approximate the detector energyresolution, decrease the sparsity of the spectral counts, andreduce the computation time needed for NMF. The time rangeof interest is 12:12:15 to 12:32:05 local time for a total of n = 2381 samples at a rate of Hz. During the collection,the aircraft made ten passes over a land-water interface at anominal altitude of ft above ground level (Figure 1).Because of the repeated passes over the land-water interface,this dataset isolates a special case of background variationoften encountered in aerial surveys. Over water and far fromland, the system essentially measures only radon and cosmicbackgrounds, plus any aircraft background, while over landit measures a combination of all background sources. Usingthis dataset, we are able to investigate the major features ofthe background and match these features to possible sourceterms.
A. NMF decomposition method
NMF reduces the dimensionality of data by solving equa-tion (1) subject to the constraint that all of the matrices arenon-negative. Essentially, NMF decomposes each spectruminto a linear combination of a pre-specified (and in ourcase, small) number of constituent spectra using non-negativecoefficients. In reality, the detection of photon events bya detector should be consistent with this formulation sincephotons themselves are non-negative and, excepting the effectsof detector dead time, detected photon events sum togetherlinearly.Interpreting the estimate ˆX = AV as the mean of thePoisson-distributed photon event counts X , the objective func-tion to minimize is the Poisson negative log likelihood Λ( X | ˆX ) = (cid:88) i (cid:88) j (cid:16) ˆ X ij − X ij log( ˆ X ij ) + log X ij ! (cid:17) . (2)There is no closed form solution to this problem so an iterativeapproach must be taken. We use the maximum likelihood update rules given in ref. [26] to perform the optimizationover both A and V , which will be referred to as training the NMF model. During the training, the rows of V areconstrained to sum to unity so that the weights A will bethe integrated number of photon counting events attributed toeach component.Two caveats to NMF are that the decomposition may not beunique and the decomposition may only find a local minimum.If there exists a k × k invertible matrix D such that V (cid:48) ≡ DV and A (cid:48) ≡ AD − are both non-negative, then A (cid:48) V (cid:48) = AV is an equivalent decomposition [25]. (Note that D itself neednot be non-negative.) A trivial solution for D is a permutationmatrix, which underscores the fact that there is no preferredordering of NMF components, unlike in PCA or SVD. Theintroduction of additional constraints on A and V in the formof regularization terms helps to break the symmetry and guidethe solutions to have desired properties [31], [34], though noconstraints are added here. Since the ordering of components isarbitrary, we chose to sort the components in decreasing orderof the variance of their weights, so that component 0 has thehighest variability, component 1 has less variability, etc. Thisordering is analogous to the ordering that results from PCAor NASVD.The choice of initialization of A and V can strongly in-fluence the particular solution that the minimization algorithmfinds, since NMF is a non-convex optimization problem andmany optimization techniques will only find the nearest localminimum [35]. Many random initializations of A and V wereinitially tried, and although some physical features (that will bediscussed later) were apparent, these components were oftennot smooth and therefore unlikely to be comparable with thespectra of different background sources.For the sake of repeatability, and also because this ini-tialization has generally led to smooth spectral components,we performed the following initialization. We made an ansatzthat NMF might at least separate the over-water backgroundfrom the terrestrial background, and to estimate the over-waterbackground we chose all spectra with total counts less than anarbitrary threshold of . times the minimum. The averagetotal counts of this set were used to initialize one columnof A , and the average unit-normalized spectrum was used toinitialize the corresponding row of V . The residual total countswere split evenly among the remaining k − components,and the corresponding rows of V were initialized to theaverage unit-normalized spectrum of the residual spectra. Thisinitialization means that the remaining k − components wouldbe initialized identically, which would result in mathematicallyidentical treatment during the optimization process. To distin-guish these k − components from each other, small positiverandom numbers less than − were added to the remaining k − components of V .As mentioned above, NMF differs from SVD in that thenumber of spectral components k must be chosen beforecalculating the decomposition. NMF models with k = 1 to components were fit to the Lake Mohave dataset, stoppingthe iterations after the change in Λ /n became less than anarbitrarily chosen threshold of − . To select the optimalnumber of components for a particular dataset, the Akaike EEE TRANSACTIONS ON NUCLEAR SCIENCE 4 A i k a k e I n f o r m a t i o n C r i t e r i o n NMF model selection
Figure 2. Akaike Information Criterion (AIC) for NMF models with differentnumbers of components for the Lake Mohave dataset. The model with threecomponents has the best (lowest) AIC. The line connecting the points hasbeen added to guide the eye.
Information Criterion (AIC) [36] was used, as in relatedwork [27]. The AIC balances the likelihood of a model withthe tendency to overfit as the number of model parametersis increased. We find that for the dataset under examinationhere, three components is the optimal number according toAIC (Figure 2).
B. NMF results
The initialized set of components for the three-componentNMF decomposition and the resulting components after NMFtraining are shown in the upper pane of Figure 3, along withthe weights in the lower pane. A comparison between theinitialized components and the final components reveals thatbesides the engineered feature of an over-water component(component 2), the two other components have divergedsignificantly in shape from each other. Some physical featuresare immediately apparent.Component 0 is characterized by prominent . and . keV lines from the decay of K and
Tl (
Thseries), respectively, and several weaker lines are apparent: . , . , . , . , and . keV, all of which arefrom the Th decay series. This evidence points to nearbyterrestrial emission as a possible source for component 0events.Component 1, on the other hand, is characterized by a highrise at the low energy end of the spectrum, peaking below keV, as well as weak line features — the only linesclearly visible are the . keV and . keV lines. Thesefeatures suggest that component 1 contains distant emissionsources with appreciable attenuation and scattering. The sharprise in the continuum below keV is also suggestiveof skyshine, the phenomenon where terrestrial emission isscattered down from the atmosphere above the detectors [22].Component 2 contains the major Rn-series lines at . , . , . , . , . , and . keV as well as thepositron annihilation line at keV. It is the dominantcomponent above keV, where a continuum of photonsdue to cosmic-ray interactions makes up the majority of thebackground [13]. Thus component 2 seems to contain radonand cosmic emission. The region around . displays anartifact from the strong K line feature in both of the othercomponents. The lower pane of Figure 3 shows the fitted weights forthe three components during the entire dataset, and Figure 4shows the same but only during the crossings of the land-water interface. The behavior of the weights at the interfacebolsters the previously mentioned spectroscopic interpretationof the NMF components. As the aircraft approaches land,component 0 rises later and more rapidly than component 1,and component 2 stays relatively constant.
C. Monte Carlo modeling of spectral components
In order to quantitatively understand the spectral shapesof the NMF components, a Monte Carlo-based model of thedetector system was used to generate spectral shapes for thepostulated sources of background. As in ref. [13], Monte Carlosimulations that can accurately model backgrounds require ageometry that extends hundreds of meters in all directions. Theauthors of ref. [13] modeled emission from the ground out toa range of m, and required an atmosphere that extended m above in order to accurately model skyshine from theground. Because our detector system is elevated above theground, we extend the ground emission out to m (2000ft) and allow air scattering up to a height and radius of m.Simulations of this scale are slow, so like ref. [13] we makesimplifying assumptions to break the simulation into smaller,more tractable parts.For each part of the simulation, a response matrix wasgenerated to relate the “input” spectrum to the “output”spectrum. Each element R ij of a response matrix R givesthe probability that an event in bin j of the input spectrumwill result in an event in bin i of the output spectrum, so theoutput spectrum is found using a matrix dot product: x out = Rx in . (3)For simplicity, the input flux distribution for all responsematrices, unless otherwise noted, is proportional to the cosineof the normal to the relevant geometry’s surface.For terrestrial emission, we generated a response matrix R terr for photon emission from within a slab of . m thicksoil by tallying the photon emission through the surface. Thesoil composition used is the U.S. average from ref. [37]. Theflux distribution that emerges is approximately proportional tothe cosine of the normal to the slab surface, as assumed in theresponse matrices and as seen in ref. [13].To simulate the transport of terrestrial photons through theatmosphere, a larger series of simulations were run using theMEGAlib Geant4-based Monte Carlo code [38], [39]. Thesesimulations involved monoenergetic gamma rays emitted froma point on a horizontal plane that represented the ground.The photons were emitted into the upward direction with anangular distribution proportional to the cosine of the anglerelative to the ground normal and were transported throughdry sea-level air until absorption or they left a cylindricalvolume of m radius and height. The tracking of photonsthat struck the ground surface were terminated. Within the air,tally volumes defined by concentric cylinders were formulatedto produce a series of track-length flux estimator tallies. Thephoton energy and azimuth and elevation angles were also EEE TRANSACTIONS ON NUCLEAR SCIENCE 5 C o un t s / k e V / s mean measurementcomponent 0component 1component 2 0 500 1000 1500 2000 2500 3000Energy (keV)10 C o un t s / k e V / s mean measurementcomponent 0component 1component 2 C o un t r a t e ( C P S ) gross countscomponent 0component 1component 2 Figure 3. NMF components and weights for a three-component decomposition of the Lake Mohave dataset. The two top panels show the componentsinitialized according to the procedure in the text (top left; components 0 and 1 are nearly identical) and the final components after optimization (top right).The bottom panel shows the gross count rate and NMF component weights. The components in the top plots are the rows of V after dividing by bin widthsand time, and then multiplying by their average weight to show their scale relative to the mean spectrum and each other. The component weights are simplythe columns of A .
300 200 100 0 100 200 300Distance from interface (m)0500010000150002000025000 C o un t r a t e ( C P S ) gross countscomponent 0component 1component 2 Figure 4. NMF component weights during the ten crossings of the land-waterinterface in the Lake Mohave dataset. Component 0 rises more rapidly thancomponent 1 when approaching land, while component 2 is approximatelyconstant. tracked in each tally volume, resulting in ◦ -spaced angularbins and 31 energy bins consisting of 30 equally spaceddown-scatter energy bins and a single full-energy bin. Theresponse due to a planar source is then inferred by leveragingsymmetry and assuming that summing tallies at different radiiis equivalent to integrating over the areal extent of the groundsource. These summed per-source photon fluxes (responses)were separately calculated at ft above a uniformly emit-ting ground plane for four geometric contributions to theoverall flux: direct responses consisting of the flux from upward-directed photons emitted within ( R dir , near ) or beyond( R dir , dist ) a “cutoff” radius r cutoff measured along the groundfrom the tally position, and skyshine responses consisting ofthe sum of all flux tally contributions from photons travelingdownward emitted within ( R sky , near ) or beyond ( R sky , dist ) thecutoff radius.For cosmic and radon backgrounds we simulated a × × m rectangular prism of atmosphere and gen-erated a response matrix for uniformly distributed isotropicemission within it and tallied photons that are incident ona × m window centered on the large face. Thisatmospheric emission matrix was denoted R atm . Radon and keV emission were propagated through this responsematrix, but the cosmic continuum was not since it is anempirical model based upon detector measurements, not thephysical production mechanism.Finally, to simulate the response of the NaI detectors toall incident background sources, a response matrix R abs wasgenerated for photons incident on the large face of a × × inch NaI detector that are absorbed by the detector.The complete response matrix for nearby emission is there-fore R near = R abs ( R dir , near + R sky , near ) R terr , (4)and the matrix for distant emission is R dist = R abs ( R dir , dist + R sky , dist ) R terr . (5) EEE TRANSACTIONS ON NUCLEAR SCIENCE 6
The complete matrix for radon and keV emission is R Rn+511 = R abs R atm , (6)while no response matrix was applied to the cosmic contin-uum. D. Results of Monte Carlo model fits
To attribute the modeled background components to theNMF components, the following analysis was performed.First, the NMF components were multiplied by their averageweights to give them an absolute scale of counts per second perkeV. This scaling was performed so that when fitting near anddistant KUT spectra to the components, not only their shapesbut also their relative magnitudes could be constrained. Thevariance of each element of this scaled V matrix was estimatedfrom the second derivatives of Λ with respect to the elementsof V (i.e., the Fisher information matrix).The nine different postulated background components werethen calculated for a given value of r cutoff and cosmic con-tinuum power-law index: • Terrestrial K (nearby and distant) • Terrestrial
U series (nearby and distant) • Terrestrial
Th series (nearby and distant) • Atmospheric
Rn series • Cosmic continuum • Cosmic keV emissionFor all sources except for the cosmic continuum, the inputspectra x in consisted of delta functions for each gamma-rayemission line weighted by branching ratio for the isotope ordecay series assuming secular equilibrium. The input spectrumfor the cosmic continuum was modeled as a power law witharbitrary normalization. All background components were thenobtained by computing the dot product with the appropriateresponse matrix from the previous section. Finally, all of thebackground components were convolved with a model of thedetector energy resolution to produce realistic line widths.In order to establish which background components con-tribute most strongly to each NMF component, a χ min-imization was performed to fit a linear combination of allnine background components simultaneously to the three NMFcomponents, for a total of 27 parameters. To find the bestcutoff radius and cosmic power-law index, the χ minimiza-tion was performed for each value of r cutoff in a grid from m to m in m increments, and for values of thepower-law index between . and . in increments of . .All coefficients were constrained to be non-negative, and thecoefficients between the nearby and distant components ofeach terrestrial type (e.g., nearby and distant terrestrial K)were constrained so that their ratio was within of theratio for uniform infinite plane emission. For all fits, the regionbelow keV was excluded because in that region materialssurrounding the detectors begin to strongly influence the shapeof the background and those materials have not been modeledhere.After performing these fits with only the constraintbetween near and distant terrestrial emission, the best-fit modelwas consistent with the basic hypothesis about the nature of the NMF components. Component 0 was dominated bythe nearby KUT background components and
Rn, whichclosely resembles the nearby
U background component.Component 1 was dominated by the distant KUT componentsand nearby
U. Finally, component 2 was dominated by
Rn, the cosmic continuum, nearby
U, and cosmogenic keV. The nearby
U and
Rn series spectra are similarin shape, which could explain the appearance of nearby
Uin the fit to NMF component 2. The best value of r cutoff was m, while the best cosmic power-law index was . .The fits were performed again but with the additionalconstraints of allowing only nearby KUT to fit component 0,distant KUT for component 1, and radon and cosmics forcomponent 2. The best overall fit was found when r cutoff =85 m and the cosmic power-law index was . . The best-fit cosmic power-law index was significantly lower than thevalue of . measured in ref. [13]. Specific details about thedifferent detectors may account for some of the difference,however such differences are unlikely to account for all ofthe discrepancy. The different elevations above sea level, aswell as the altitude above ground level, may also affect thepower law observed. Differences in how radon was modeledcould also contribute. The nature of this discrepancy remainsunknown.Figure 5 shows the results of these constrained fits. Theskyshine portion from all terrestrial components was summedtogether and plotted separately to show the portion of thecomponent that is estimated to come from skyshine. Forcomponent 0, the skyshine contribution makes up a maximumof about of the spectrum below 100 keV, while forcomponent 1 skyshine makes up two-thirds of events below keV and more than half of all events below keV. Com-ponent 2 is dominated by radon, with smaller contributionsfrom the cosmic continuum and keV emission. The sum ofall the weighted NMF components, which is approximately themean measured spectrum, is also shown in Figure 5, togetherwith the sums of all the fitted background components.The NMF components and the background model fits are ingeneral agreement and provide strong evidence for the iden-tification of component 0 with nearby emission, component 1with distant emission, and component 2 with radon and cosmicemission. E. Comparison to NASVD
Finally, we compare the NMF components to componentsderived using NASVD. To perform NASVD, one scales eachelement of X to have unit variance and then performs singularvalue decomposition (SVD) on the resulting matrix [21]. Toestimate the mean of each bin, NASVD assumes that theaverage spectral shape is approximately the same. Defining s to be the average of the rows of X that is then normalized tosum to (i.e., the average unit-normalized spectrum): s j ≡ (cid:80) i X ij (cid:80) i (cid:80) j X ij , (7) EEE TRANSACTIONS ON NUCLEAR SCIENCE 7
Energy (keV)10 C o un t s / s / k e V NMF component 0simulation totalskyshine totalK-40 (near)U-238 series (near)Th-232 series (near) F i t r e s i d u a l () Energy (keV)10 C o un t s / s / k e V NMF component 1simulation totalskyshine totalK-40 (dist)U-238 series (dist)Th-232 series (dist) F i t r e s i d u a l () Energy (keV)10 C o un t s / s / k e V NMF component 2simulation totalRn-222 seriescosmic continuumcosmogenic 511 F i t r e s i d u a l () Energy (keV)10 C o un t s / s / k e V NMF totalsimulation totalskyshine totalK-40 (near)K-40 (dist)U-238 series (near)U-238 series (dist)Th-232 series (near)Th-232 series (dist)Rn-222 seriescosmic continuumcosmogenic 511 F i t r e s i d u a l () Figure 5. Monte Carlo background components fit to the three NMF components (top left, top right, and bottom left) and the sum of all three components(bottom right) from the Lake Mohave dataset. These fits provide strong evidence for the identification of component 0 with nearby emission, component 1with distant emission, and component 2 with radon and cosmic emission. The skyshine contributions from the terrestrial sources are shown separately. and also defining the sum of the counts in each spectrum (i.e.,the gross counts) to be c : c i ≡ (cid:88) j X ij , (8)NASVD uses c i s j as an approximation to the mean of mea-surement X ij and forms the unit-variance matrix X (cid:48) : X (cid:48) ≡ diag( c ) − / X diag( s ) − / . (9)Then SVD is used to decompose X (cid:48) : X (cid:48) = UΣW T , (10)where U is an n × n unitary matrix, Σ is an n × d diagonalmatrix of the singular values, and W is a d × d unitary matrix.To remove noise, the highest k singular values are chosen andthe rest are culled to form the n × k matrix Σ k and the d × k matrix W k . Now defining A ≡ diag( c ) / UΣ k (11) V ≡ W Tk diag( s ) / (12)we obtain the NASVD solution to equation (1). A consequenceof NASVD is that the first column of UΣ k is √ c and the firstcolumn of W k is √ s , so therefore the first column of A is c and the first row of V is s [40]. The remaining componentsare additive perturbations to the mean spectrum in order ofdecreasing variance, and the row-wise sums of V beyond thefirst component are zero. (This latter fact is a result of theorthogonality of W k . Letting w j be the j th column of W k , then since the first column of W k is w = √ s , the sum ofthe j th row of V is w j · √ s = w j · w = δ j .)Figure 6 shows the results of an NASVD decompositionfor the Lake Mohave dataset. Components 0–2 account forover of the variance, while components 3 and higherhave much smaller variances than the first three and do notshow coherent spectral features. Components 1 and 2 clearlydisplay features that can be associated with spectroscopicphenomena, such as photopeaks, but because they take onnegative values, each component cannot be the direct resultof background source emission. For example, increasing theweight of component 1 decreases the relative contribution ofradon and cosmics while increasing the contribution of nearby K and
Tl emission. Likewise, increasing the weight ofcomponent 2 increases the contribution from distant emissionrelative to nearby emission. As a further comparison to theNMF decomposition, Figure 7 shows the weights of the firstthree NASVD components during the ten crossings of theland-water interface and is the equivalent of Figure 4. At theinterface, the weights for components 1 and 2 change sign toaccount for the rapidly changing shape of the spectrum in thatregion.In order to quantitatively compare the abilities of NASVDand NMF to represent the same spectra, we first note thatthe gross counts of both NASVD- and NMF-reconstructedspectra exactly match the gross counts of the measured spectra.For NASVD, this property follows from the fact that the firstcolumn of A are the gross counts c , and the sum of the j throw of V is δ j . For NMF, this property follows from the EEE TRANSACTIONS ON NUCLEAR SCIENCE 8 C o un t s / k e V / s component 0component 1component 2component 3component 4 C o un t r a t e ( C P S ) component 0 × 0.1component 1component 2component 3component 4 Figure 6. NASVD components for the Lake Mohave data (top) and component weights (bottom). Spectral components have been multiplied by their maximumweight to show their scale relative to each other. In the bottom plot, the component 0 weights have been scaled by a factor of . for ease of comparison tothe others. The first NASVD component is always the average spectral shape and its weights are the gross counts, while the following components representvariations from the average shape. Components 0–2 account for over of the variance of the dataset while components 3 and higher largely containstatistical noise.
300 200 100 0 100 200 300Distance from interface (m)1000010002000 C o un t r a t e ( C P S ) component 0 × 0.1component 1component 2 Figure 7. Weights of the first three NASVD components during the tencrossings of the land-water interface in the Lake Mohave dataset (cf. the NMFweights in Figure 4). The weights for component 0 are the gross counts, whilecomponents 1 and 2 change the relative importance of radon and cosmics anddistant emission, respectively. The component 0 weights have been scaled bya factor of . for ease of comparison to the others. fact that when maximizing Poisson likelihood, any model thathas an effective scale factor will be fit so that the sum of themodel matches the measured gross counts.Since both methods exactly reconstruct each spectrum’sgross counts, it is also important to compare how well theyfit the shape of the spectrum. For this purpose we introducethe Spectral Anomaly Detection (SAD) metric, the L2 normbetween the i th measured spectrum and the i th reconstructed spectrum [41], [27]: SAD i ≡ (cid:80) j || X ij − ˆ X ij || (cid:80) j || X ij || , (13)where the denominator normalizes the metric to a mean ofunity. We evaluated SAD for NMF and NASVD reconstruc-tions with three components and histogrammed the resultsin Figure 8. The distribution of the difference in the SADmetric between the NASVD and NMF reconstructions ismuch narrower than the distribution of either SAD metric,indicating that the two methods reconstruct the measuredspectra with similar fidelity. This result is not surprising sinceboth models have the same linear form with similar numbers offree parameters, though NMF has additional (non-negativity)constraints. Notably, the histogram of SAD metric differenceshas an average value greater than zero, which indicates thatNMF yields slightly better reconstructions than NASVD, atleast according to the SAD metric. The average value ofthe SAD metric difference becomes negative for this dataset(i.e., favoring NASVD) only once the number of NASVDcomponents used in the reconstruction has been increased tonine.III. A PPLICATION TO S AN F RANCISCO B AY A REASURVEYS
To demonstrate the application of NMF decompositionunder more general conditions, two survey datasets fromthe San Francisco Bay Area (henceforth, Bay Area) were
EEE TRANSACTIONS ON NUCLEAR SCIENCE 9 C o un t s p e r b i n SAD
NASVD
NMF
NASVD
SAD
NMF
Figure 8. Spectral Anomaly Detection (SAD) metric for spectra reconstructedusing NASVD and NMF with 3 components. The red histogram shows thedifference in the metric between the two methods. This difference is small,indicating that they reconstruct the measured spectra with similar fidelity. analyzed using the techniques of the previous section. Thesedatasets were part of a larger survey campaign performedin August 2012 [42]. The first set was taken at the Port ofOakland on 30 August 2012 and consisted of 5,687 -Hzspectra, and the second dataset was collected in Pacifica on28 August 2012 and consisted of 4,505 -Hz spectra. The datawere taken with the same NaI(Tl) detector system as the LakeMohave data, and the same energy range and binning schemewere used. Both the data rate ( Hz) and the nominal flightaltitude above ground level ( ft) differed from the LakeMohave dataset ( Hz and ft, respectively). Both surveyscontain land-water interfaces like the Lake Mohave survey,but both also contain roads and man-made structures. ThePacifica survey is also of special interest because it containsrugged terrain and some geological features with substantiallylower potassium content than typical (and therefore weak Kemission), providing additional spectral variability.
A. Anomaly detection and removal
For the Bay Area datasets, the Spectral Comparison RatioAnomaly Detection (SCRAD) method [43], [44], [7] was ap-plied to the data to identify and remove any spectral anomaliesthat could affect the NMF decomposition. To use SCRAD oneneeds to bin each spectrum into a series of predeterminedbins, track estimates of the means and covariances of the binsover time using an exponentially weighted moving average(EWMA), and then construct a vector of spectral comparisonratios (SCRs) from the binned counts. SCRAD relies on thebackground estimate being accurate (we used an EWMA witha parameter of . ) and each bin being large enough that itsstatistics are Gaussian, not Poisson.For SCRAD, 13 non-overlapping spectral windows wereused with boundaries of , , , , , , , , , , , , , and keV. These windowswere chosen to span the entire energy range and so that nocounts in any bin were fewer than in order to ensure theGaussian approximation was statistically valid. No effort wasmade to optimize the windows for any particular types ofanomalies, nor was nuisance rejection implemented (i.e., N- SCRAD [45]). The SCR distance metric D SCR was used asthe test statistic [44], and its expected chi-squared distributionwith degrees of freedom was used to set a threshold of σ significance.Using this metric, three anomalies were found in the Portof Oakland dataset. These anomalies, at indices , ,and , can be seen in Figure 10 where they are markedwith red arrows. The anomaly at index lasted at least seconds and is of unknown origin. It consisted of a hardcontinuum extending beyond MeV. The anomaly at index consisted of an increase in counts around keV andoccurred when the aircraft was near the position of the firstanomaly on the following flight line, so it is likely down-scattered photons from the first anomaly. The third anomalyat index was likely a
Tc source due to elevated countsup to keV. The data associated with these anomalies wereexcluded from the NMF training and are not shown on themaps in Figure 11.
B. Results of NMF decompositions
NMF was initialized and trained on the Bay Area datasetsin the same manner as the Lake Mohave dataset. NMF modelswith k = 1 to components were trained, and AIC wasagain used to determine the optimal number of components todescribe the dataset, resulting in k = 3 components for boththe Port of Oakland and Pacifica datasets. Figure 9 shows theresulting NMF components for each dataset, and Figure 10shows the weights for the Port of Oakland dataset.NMF decompositions are not expected to be identical in allsituations since they are fitted to the training data provided.For example, the Pacifica survey was over rugged terrain withregions of abnormally low K, while the Port of Oaklandsurvey was over flat terrain but with man-made materials,which could give rise to sharp changes in KUT concentrations.The datasets may also have different ambient backgrounds,which would affect any terms that are approximately constant.In general both NMF decompositions are qualitatively sim-ilar to each other and to the Lake Mohave decomposition(Figure 3), aside from some notable differences. Component 0on average makes up a smaller fraction of the total spectrum,which may be due to greater attenuation of nearby emission atthe higher altitude. One subtle difference in component 0 fromPacifica is that there is an indication of a weak
Cs peak at keV from 1950s weapons-testing fallout, perhaps due toless surface disturbance than the Port of Oakland and greaterdeposition than Lake Mohave due to higher yearly rainfall. Incomponent 1, the K peak at . keV is less prominent(Oakland) or not at all present (Pacifica). Component 2 inboth Bay Area datasets now includes a clear K peak, whichcould be due to different aircraft background and the presenceof potassium in the salt water of the Pacific Ocean and theSan Francisco Bay.Maps were generated for both datasets by plotting the NMFcomponent weights on a map without any further correc-tions (Figure 11). As expected from the reasoning developedin Section II, component 0 has more well defined land-waterboundaries than component 1, consistent with component 0
EEE TRANSACTIONS ON NUCLEAR SCIENCE 10 C o un t s / k e V / s mean measurementcomponent 0component 1component 2 0 500 1000 1500 2000 2500 3000Energy (keV)10 C o un t s / k e V / s mean measurementcomponent 0component 1component 2 Figure 9. NMF components for a three-component factorization fitted to the Port of Oakland data (left) and Pacifica (right). These decompositions have astrong resemblance with each other and with the Lake Mohave decomposition (Figure 3). A red arrow marks the presence of the keV line from
Csin component 0 from Pacifica. C o un t r a t e ( C P S ) gross countscomponent 0component 1component 2 D S C R t e s t s t a t i s t i c Figure 10. The gross count rate and NMF component weights for the Port of Oakland dataset (top). The bottom plot shows the SCRAD test statistic usedto identify the anomalies in the dataset, and the red line indicates a threshold of σ significance. The red arrows mark the locations of the anomalies. representing nearby terrestrial emission and component 1representing distant terrestrial emission. Component 2 is flatfor the entirety of both surveys, which is consistent withcosmic and radon emission. C. Modeling the source terms
The NMF components were modeled in a similar way asthe Lake Mohave dataset (Section II-C), with the primary dif-ference being that the response matrices R dir , near , R sky , near , R dir , dist , and R sky , dist were adjusted for the higher altitudeusing the same suite of simulations. In an attempt to modela potential K background from the aircraft, an extra Ksource was fit to component 2 that consisted of R abs appliedto a K spectrum. The best value of r cutoff was determinedto be m for the Port of Oakland and m for Pacifica (cf. r cutoff = 85 m for the Lake Mohave dataset at ft).The optimal power-law index was found to be . for bothdatasets, significantly larger than the Lake Mohave value( . ) but closer to the measurement in ref. [13] ( . ). It isunknown why there is such a large discrepancy between theLake Mohave power-law index and the Bay Area power-lawindex, but it could point to a deficiency in the modeling of theradon, cosmics, or both.The results of fits between the Monte Carlo backgroundcomponents and NMF components for the Port of Oakland areshown in Figure 12. These fits were constrained in the sameway as the fits to the Lake Mohave components. Once again,the modeled backgrounds are a qualitatively good fit to theNMF components, providing evidence for the identification ofcomponent 0 with nearby emission, component 1 with distant EEE TRANSACTIONS ON NUCLEAR SCIENCE 11
Figure 11. Maps of the weights for each of the three NMF components in the Port of Oakland (left) and Pacifica (right) datasets. For both NMF decompositions,component 0 (top) has sharper land-water boundaries than component 1 (middle), while component 2 (bottom) remains relatively constant. (Map imagery:Google.) emission, and component 2 with radon and cosmic emissionat this altitude.As in the Lake Mohave fits, the model for component 0 hasan excess continuum below keV, although these energiescould easily be attenuated by a small amount of material nearthe detectors. Unlike the Lake Mohave fits, the model forcomponent 2 also has an excess below keV due to thelarger contribution of the cosmic continuum relative to radon(possibly due to the higher altitude above ground level). Thisexcess could also be attenuated by material near the detectors,which is not included in this model.IV. D
ISCUSSION
In this work we have presented the application of NMF tothe decomposition of gamma-ray spectra from aerial surveys,with a particular focus on identifying background source termsfor the resulting spectral components. Since NMF uses the dataitself to derive the primary spectral components, it has theadvantage over FSA of being able to adapt to spectral shapes that have not been included in modeled components. Whencompared to NASVD, NMF is able to reconstruct spectra toa similar fidelity while preserving both non-negativity andmaximizing Poisson likelihood, and, in at least the casespresented here, the components appear to have a plausiblephysical origin.Here we have focused on aerial data that includes water, afeature that isolates the radon and cosmic background com-ponents and favors a three-component NMF decomposition.If the data do not include time over water, it may be harderto separate radon and cosmics from other background. Forthese situations, other methods such as altitude spirals, wherethe aircraft flies repeated patterns in the same location atincreasing altitudes, may help separate cosmics and radonfrom other background. Also, a validated Monte Carlo modelcould potentially be used to initialize an NMF decompositionin order to guide it toward physically meaningful components.A feature of NMF decompositions in the cases presentedhere is that NMF can approximately separate distant terrestrial
EEE TRANSACTIONS ON NUCLEAR SCIENCE 12
Energy (keV)10 C o un t s / s / k e V NMF component 0simulation totalskyshine totalK-40 (near)U-238 series (near)Th-232 series (near) F i t r e s i d u a l () Energy (keV)10 C o un t s / s / k e V NMF component 1simulation totalskyshine totalK-40 (dist)U-238 series (dist)Th-232 series (dist) F i t r e s i d u a l () Energy (keV)10 C o un t s / s / k e V NMF component 2simulation totalRn-222 seriescosmic continuumcosmogenic 511K-40 (aircraft) F i t r e s i d u a l () Energy (keV)10 C o un t s / s / k e V NMF totalsimulation totalskyshine totalK-40 (near)K-40 (dist)U-238 series (near)U-238 series (dist)Th-232 series (near)Th-232 series (dist)Rn-222 seriescosmic continuumcosmogenic 511K-40 (aircraft) F i t r e s i d u a l () Figure 12. Monte Carlo background components fit to the three NMF components (top left, top right, and bottom left) and the sum of all three components(bottom right) from the Port of Oakland dataset. The skyshine contributions from the terrestrial sources are shown separately. emission from other background sources using the data alone.Evidence for this separation has also been seen in vehicle-borne gamma-ray data [32]. By exploiting this separation,NMF can potentially improve the resolution of aerial surveymaps by disentangling distant emission from nearby emission.In addition, by finding background components with phys-ical origins, NMF could be leveraged in new algorithms foranomaly detection and background estimation. NMF has al-ready been shown to be competitive with PCA-based methodsfor spectral anomaly detection, which may be due to its accu-rate treatment of Poisson statistics and ability to be consistentwith physics [27]. The different temporal variability of thecomponents could be exploited by Kalman filters or low-passfilters to find anomalous behavior. Even in the data presentedhere from the Port of Oakland (Figure 10), the anomaliesfound by the SCRAD metric can be easily identified as briefspikes in the NMF component 1 count rate. Since component 1is believed to arise from distant emission sources, it shouldnot exhibit such high frequency behavior. Component 2 alsoincreases during the large anomaly at index , which is anunusual departure from its relatively constant count rate.With the ability to decompose the low-energy region of thespectrum where skyshine and down-scatter are most signifi-cant, NMF could be used by algorithms intended to identifyanomalies below keV. This regime is especially of interestwhen searching for spectral anomalies due to nuclear material,since many important isotopes have gamma-ray lines below keV (e.g., [46]), however that region of the spectrum isnotoriously difficult to model.One drawback with using NMF is that the training can
TABLE IC
OMPUTATION TIMES FOR
NMF
ON THE L AKE M OHAVE DATASET . Method Components Training Single spectrumNASVD N/A . s N/ANMF k = 1 0 . s . msNMF k = 2 65 . s . msNMF k = 3 107 . s . msNMF k = 4 389 . s . msNMF k = 5 1 , . s . msbe computationally intensive — the k = 3 model tookabout s running on four . GHz cores, however thematrix computations can easily be split across any number ofavailable cores. TABLE I shows the duration of the trainingtime for different NMF models trained on the Lake Mohavedataset, including the time it takes to calculate the weightsfor a single spectrum after the training has been completed.The time to calculate the NASVD decomposition is shown forcomparison. The length of training time is set largely by thearbitrary convergence criterion used here, and NMF models ingeneral may not need to meet the same criterion to be useful.Even though the training can take a significant length oftime, once the NMF model has been trained for a givendetector and environment and V is fixed, the subsequentcalculation of weights for previously unseen spectra is fast.This application highlights another distinction between FSA,NMF, and NASVD — typically, FSA is used to decomposespectra that have not yet been measured, while NASVD isused to decompose only previously measured spectra (although EEE TRANSACTIONS ON NUCLEAR SCIENCE 13 see ref. [24] for a real-time application of NASVD). NMFhas been presented here retrospectively like NASVD, but bothNMF and NASVD components can be fit to future measure-ments under the assumption that the components derived fromthe training set are also effective at reducing the dimensionalityof future data. Using a trained NMF model to decomposespectra has been shown to be an effective tool for spectralanomaly detection in a mobile detection system [27].Finally, NMF is a potentially powerful framework for anyanalysis involving variable gamma-ray backgrounds since itprovides a consistent description for both the mathematicaland physical characteristics of the measurements. As such,NMF provides a natural framework for data fusion, such asincorporating non-radiological contextual information aboutthe environment. Some work on vehicle backgrounds hasalready shown evidence for the ability to connect NMF com-ponents to non-radiological contextual information (i.e., videoimages) [32]. As a simple example of data fusion for aerialmeasurements, having knowledge of the distance to a land-water interface would provide useful information about thebehavior of the distant and nearby background components.This area has much promise for further research.V. A
CKNOWLEDGMENTS
This work has been supported by the US Department ofHomeland Security, Domestic Nuclear Detection Office, underIAA HSHQDC-11-X-00380. This support does not constitutean express or implied endorsement on the part of the Govern-ment.This work was also performed under the auspices of theUS Department of Energy by Lawrence Berkeley NationalLaboratory under Contract DE-AC02-05CH11231. The projectwas funded by the US Department of Energy, National NuclearSecurity Administration, Office of Defense Nuclear Nonpro-liferation Research and Development (DNN R&D).We would like to thank Ren Cooper and Joseph Curtis fortheir helpful comments and Erika Suzuki for proofreading thismanuscript. R
EFERENCES[1] R. Moxham, “Airborne radioactivity surveys in geologic exploration,”
Geophysics , vol. 25, no. 2, pp. 408–432, Apr. 1960. [Online]. Available:https://library.seg.org/doi/abs/10.1190/1.1438711[2] J. H. Galbraith and D. F. Saunders, “Rock classification by charac-teristics of aerial gamma-ray measurements,”
Journal of GeochemicalExploration
Journal of Environmental Radioactivity
Geophysics , vol. 40, no. 3, pp. 503–519, Jun. 1975.[Online]. Available: http://library.seg.org/doi/abs/10.1190/1.1440542[6] R. Grasty, J. Glynn, and J. Grant, “The analysis of multichannelairborne gammaray spectra,”
Geophysics , vol. 50, no. 12, pp. 2611–2620, Dec. 1985. [Online]. Available: https://library.seg.org/doi/abs/10.1190/1.1441886 [7] R. S. Detwiler, D. M. Pfund, M. J. Myjak, J. A. Kulisek, andC. E. Seifert, “Spectral anomaly methods for aerial detection usingKUT nuisance rejection,”
Nuclear Instruments and Methods in PhysicsResearch Section A: Accelerators, Spectrometers, Detectors and Associ-ated Equipment
Geophysics , vol. 47, no. 1, pp. 117–126, Jan. 1982.[Online]. Available: https://library.seg.org/doi/abs/10.1190/1.1441273[9] B. Minty, “Airborne gammaray spectrometric background estimationusing full spectrum analysis,”
Geophysics , vol. 57, no. 2, pp. 279–287,Feb. 1992. [Online]. Available: http://library.seg.org/doi/abs/10.1190/1.1443241[10] B. Minty, P. McFadden, and B. Kennett, “Multichannel processingfor airborne gammaray spectrometry,”
Geophysics , vol. 63, no. 6, pp.1971–1985, Nov. 1998. [Online]. Available: https://library.seg.org/doi/abs/10.1190/1.1444491[11] H. K. Aage, U. Korsbech, K. Bargholz, and J. Hovgaard, “Anew technique for processing airborne gamma ray spectrometrydata for mapping low level contaminations,”
Applied Radiation andIsotopes
AGSOJournal of Australian Geology & Geophysics , vol. 17, no. 2, pp. 39–50,1997.[13] G. A. Sandness, J. E. Schweppe, W. K. Hensley, J. D. Borgardt,and A. L. Mitchell, “Accurate Modeling of the Terrestrial Gamma-Ray Background for Homeland Security Applications,” in
Nuclear Science andEngineering
Journal of Environmental Radioactivity
Science of The Total Environment
Canadian Journal of Earth Sciences
Proceedings from the Sixth Topical Meeting on Emergency Preparedness& Response , 1997.[22] J. Hovgaard and R. Grasty, “Reducing Statistical Noise in AirborneGamma-Ray Data Through Spectral Component Analysis,” in
Proceed-ings of Exploration 97 , 1997, pp. 753–764.[23] B. Minty and P. McFadden, “Improved NASVD smoothing of airbornegamma-ray spectra,”
Exploration Geophysics , vol. 29, no. 3/4, pp.516–523, Dec. 1998. [Online]. Available: http://library.seg.org/doi/abs/10.1071/EG998516[24] J. A. Kulisek, J. E. Schweppe, S. C. Stave, B. E. Bernacki, D. V.Jordan, T. N. Stewart, C. E. Seifert, and W. J. Kernan, “Real-timeairborne gamma-ray background estimation using NASVD with MLEand radiation transport for calibration,”
Nuclear Instruments andMethods in Physics Research Section A: Accelerators, Spectrometers,Detectors and Associated Equipment , vol. 784, pp. 287–292, Jun.
EEE TRANSACTIONS ON NUCLEAR SCIENCE 14
Environmetrics , vol. 5, no. 2, pp. 111–126, Jun. 1994.[Online]. Available: https://onlinelibrary.wiley.com/doi/abs/10.1002/env.3170050203[26] D. D. Lee and H. S. Seung, “Learning the parts of objects bynon-negative matrix factorization,”
Nature
IEEE Transactions on Nuclear Science
IEEE Transactions on Medical Imaging , vol. 23, no. 12, pp.1453–1465, Dec. 2004.[30] C.-Y. Liou and K.-D. Ou Yang, “Unsupervised Classification of Re-mote Sensing Imagery With Non-negative Matrix Factorization,” in
InICONIP , 2005.[31] V. P. Pauca, J. Piper, and R. J. Plemmons, “Nonnegative matrixfactorization for spectral data analysis,”
Linear Algebra and itsApplications
Nuclear Instruments andMethods in Physics Research Section A: Accelerators, Spectrometers,Detectors and Associated Equipment
MIPPR 2007: Remote Sensing and GIS Data Processingand Applications; and Innovative Multispectral Technology andApplications
Advances in Neural Information Processing Systems13 , T. K. Leen, T. G. Dietterich, and V. Tresp, Eds. MIT Press,2001, pp. 556–562. [Online]. Available: http://papers.nips.cc/paper/1861-algorithms-for-non-negative-matrix-factorization.pdf[36] H. Akaike, “A new look at the statistical model identification,”
IEEETransactions on Automatic Control
New AstronomyReviews
Nuclear Instruments and Methods in Physics ResearchSection A: Accelerators, Spectrometers, Detectors and AssociatedEquipment
Airborne gamma-ray spectrometry . Kgs. Lyngby, Den-mark: Technical University of Denmark (DTU), Mar. 1998.[41] K. Miller and A. Dubrawski, “Gamma-Ray Source Detection With SmallSensors,”
IEEE Transactions on Nuclear Science
IEEE Transactions on Nuclear Science , vol. 54, no. 4, pp.1232–1238, Aug. 2007.[44] K. D. Jarman, R. C. Runkle, K. K. Anderson, and D. M. Pfund,“A comparison of simple algorithms for gamma-ray spectrometersin radioactive source search applications,”
Applied Radiation andIsotopes
Applied Radiationand Isotopes