This chapter presents some recent data processing developments associated with radiation monitoring systems. Radiation monitors have to continuously provide count rate estimations with accuracy and precision. A filtering technique based on a Centered Significance Test coupled with a Brown’s double exponential filter has been developed and used in compensation measurement and moving sources detection schemes.
- data processing
- nuclear counting
- radiation monitor
- signal processing
- frequentist inference
During the last decades, ionizing ray detectors have grown in performance, thanks to digital electronics developments (ADC and FPGA), allowing for an advanced processing of nuclear impulse signals. It is also noteworthy that this field has favored the development of real-time processing algorithms dealing with count rate data.
The architecture of a typical nuclear measurement system is presented in Figure 1. It can be divided in four parts:
The detector part contains the physical sensor (noble gas, scintillation material, and semiconductor) in which radiation interacts with matter. A conversion unit (preamplifier or photo-converter) converts the induced charges or photons in amplified voltage pulses. In the case of gas or semiconductor-based sensors, a high voltage is required to polarize the medium, and a low voltage is needed to supply active components of the preamplifier. In the case of scintillators, a high voltage supplies the photomultiplier.
Front-end electronics is composed of an analog filter, an analog-to-digital converter (ADC), and a digital filter. An analog shaping filter can be used to adapt the signal before digital conversion (dynamic range and respect to Shannon rules), and/or to maximize signal-to-noise ratio (SNR). The ADC digitalizes the signal with a given frequency and resolution.1 This digital signal is processed into a fast electronic component, typically a microcontroller, or a field-programmable gate array (FPGA). The embedded firmware has to comply with the very high-frequency of the ADC output with a processing period in the range of 1–10 ns. The algorithms implemented in the firmware perform the pulse processing, which mainly consists in triggering, first digital filtering for SNR maximization, stabilizing the baseline, estimating the dead time, and counting a number of pulse events N over a period of time ∆τ. This general description is not exhaustive and a variety of architectures is conceivable, depending on the mix between analog and digital processing. Though modern trends tend to favor digital filtering, analog filtering can still be retained to comply with cost reduction or embedded strategies. For instance, the front-end associated with a scintillator can directly digitalize the output voltage of a photomultiplier using a 500 MS/s ADC and process the pulses using a FPGA (notably when pulse shape discrimination is needed). On the other hand, the initial signal can be filtered using analog components (trapezoidal filtering), before digitalization with a 10 MS/s ADC and count processing with a microcontroller.
An interface is built on a computer connected with the front-end electronic card. The software reads, at each given time interval ∆t, a new count value Ni according to a defined communication protocol. This second processing can be divided in two parts: filtering of the count rate signal and displaying. The period ∆t has to be chosen in compliance with continuous measurement requirement, typically close to retinal persistence of 0.1 s. We can highlight here, the quantitative difference between the digital pulse processing coded in VHDL or Verilog for the FPGA (very fast), and the count rate processing which can be coded in C/C++ into a microcontroller or/and the PC interface (7–8 orders of magnitude slower). In compact systems, the count rate processing is usually incorporated into the firmware while, in larger systems, the count rate processing is remotely implemented in the PC.
This chapter will not address pulse processing techniques, for which details can be found in [1–3], but presents some recently developed techniques to process count rate signal using frequentist inference. Bayesian inference can also be implemented to process count rate as for instance for gamma spectrum unfolding or photon-limited imaging filtering [4, 5]. These are very efficient to accurately processed nuclear counting data, but become unsuited of online applications. After describing the theoretical model of the counting process, a smoothing technique will be presented as a fundamental building block, ensuring an online and adaptive filtering of the signal. The issue of composite measurements will then be addressed with a method allowing improving metrological reliability for particle discrimination (compensation technique). Finally, the use of detectors in a network to address moving source detection will be developed.
2. Nuclear counting model
Nuclear disintegration can occur following different processes depending on the A/Z ratio of the concerned isotope. Major disintegration processes read: β−, β+, ε, α, and spontaneous fission decay, presented in the following nuclear equations, where X is the mother nucleus and Y the daughter nucleus
Subsequently, the daughter nucleus is, most of the time, released in an exited state and usually reaches its fundamental level by gamma-ray emission:
According to the detector type, β−, β+, α, n, or γ particles are detected and counted. The required time τd for an unstable nucleus to decay is undetermined, and takes its value in an exponential distribution whatever the time lap between its creation and the observation is (memoryless phenomenon). The probability distribution p(τd = t′) of the decay instant, where t′ and λ are, respectively, the observation instant and the decay constant of the nucleus, is given by:
The observation of an unstable nucleus over a time t forms a Bernoulli trial in which two results can be observed: the nucleus has decayed or the nucleus has not decayed. A probability of p and 1 − p can be, respectively, associated with each branch of the trial for every instant as illustrated in Figure 2.
The probability pX → Y(t) to observe a disintegration of the mother nucleus X → Y before time t (t′ = 0 being the start of the observation) is obtained as:
In a radioactive source containing a population of NX unstable nuclei, the decay of an individual nucleus does not impact the decay of the others. The Bernoulli trial is therefore repeated NX times (Figure 3) during the observation time t, and the number of observed decays n is described by a Binomial law such as:
In practice, NX is very large and pX → Y(t) is usually very small (1/λ ≫ t). In these conditions, the Binomial law converges toward a Poisson law such as:
Expectation and variance of the number of decays are equal to NXλt. The number of counts N measured before observation time t is obtained by weighting Eq. (9) with the detection efficiency ε and the probability η of the detected particle to be emitted during decay. The expected count rate ρ = εηNXλ thus becomes the parameter of the distribution of measured count values before t:
At each time ti, sampled such as t0 = 0 & ∀ i ≥ 1, ti = ti − 1 + ∆t, the raw estimation of the count rate ρ(ti) = ρi is provided by measuring Ni, which is a time-dependent random variable taking its values in a Poisson distribution such as:
A challenge in radiation monitoring is to provide count rate estimation ρi at each time ti maximizing both precision σ(ρi) → 0 and accuracy σ(ti) → 0. Algorithmic techniques to meet this expectation are discussed in the next section.
3. Count rate smoothing
The aim of smoothing algorithms is to improve the estimation of ρi, originally defined as: . This improvement can be achieved by using past values Ni − 1, Ni − 2, … recorded in a memory according to the assumption that p(ρi|Ni, Ni − 1, Ni − 2, …) is more precise than p(ρi|Ni). If we consider, in a first approach, a constant count rate ρ, the estimator which maximizes the likelihood of a homogenous Poisson process is the average :
where m + 1 is the temporal depth of the filter and ϑ is a kernel function in which each ϑj, 1 ≤ j ≤ i equal to one. According to the property of equality between variance and expectation, the associated variance can be estimated as:
The relative stochastic uncertainty is inversely proportional to the square root of the historical depth m + 1.
In practice, counting processes are not homogenous (ρ is not constant). In this case, it is important to provide an estimate of the time effectively corresponding to the current count rate estimate . Because the sampling times ti over the temporal depth m + 1 are identically weighted, the estimate from Eq. (12) is associated to a time estimate with a temporal precision . We therefore see that can only be minimized to the detriment of , leading to a degradation of accuracy when ρ is varying. One way to address this issue is to actualize the temporal depth mi after every count rate estimation. The optimal value for mi is a function of the temporal behavior of the rate ρ at the time or .
First approaches consist in the implementation of preset count filter providing a fixed variance , or finite impulse response (FIR) filters in which a kernel function ϑ is used to assign more weight to recent than older count values such as Eq. (12) become:
The algorithm translation of the actualization of mi is the building of infinite impulse response (IIR) dedicated to nuclear counting . Such nonlinear filtering requires a hypothesis test to detect the changes in count rate ρ. The null hypothesis H0 and the detection hypothesis H1 are formalized as follows:
In a first approach , a sequential probability ratio test (SPRT) has been assessed under the assumption that θ1 is a known value. Later, generalized tests in which θ1 is an unknown parameter have been introduced, notably the generalized likelihood ratio test (GLR)  and the centered significance test (CST) . In these change detection algorithms, several estimations of the current count rate are calculated using different temporal depths k such as:
In the rest of the discussion, we will conventionally use notation to designate both the underlying random variable and its actual values.
In the CST test, the vector is scanned to find a potential change in the true rate ρ. For every temporal depth k ∈ ⟦1; mi⟧, the difference between count rate estimations is the quantity which will be tested for significance.
The method is based on a comparison between actual and expected distributions of under H0 and H1, respectively . The distribution of is the difference between two weighted Poisson distributions ℘
The expectation between times (i − mi)∆t and i∆t. Moreover, we will make use of assumption , as common in nuclear counting experiments will finite statistics . The variances associated with both uncorrelated random processes are summed to obtain a cumulative standard deviation for . According to the equality between expectation and variance, we obtain:
We will note , the distribution of the difference random variable with its first and second order moments.
Under H0, (cf. left curve in Figure 4). A decision threshold (DT) is determined in compliance with a given risk of false detection . DT is defined in the following formula, where is the quantile of the error function (err) with a confidence level:
In practice, for embedded implementations, it is impossible to sample and interpolate distributions for every value of i and k. Moreover, when and are large enough, distribution may be approximated as , with the Normal law. Under this assumption, for every value of i and k, and , where Q1 − α is a quantile of and err becomes:
where Φ is the cumulative distribution function of the centered Normal law.
As illustrated in Figure 4, DT can be calculated thanks to the weighting of the standard deviation by Q1 − α such as:
If , hypothesis H0 is accepted with a confidence level equal to 1 − α.
Under H1, (cf. right curve in Figure 4), where DL is defined as the detection limit. DL is determined in compliance with a given risk of non-detection β = p(H0| H1). DL is obtained in the following formula:
As illustrated in Figure 4, DL can be calculated thanks to the weighting of the associated standard deviation by the quantile Q1 − β of the error function erf such as:
An equivalent confidence level 1 − α = 1 − β = γ is considered, and Eq. (23) is solved recursively such as:
When y → ∞,
If , hypothesis H1 is accepted with a confidence level 1 − β = γ.
The number Li of significant changes recorded into memory is used to calculate the next value of temporal depth mi + 1:
If Li = 0, true rate ρ is considered to remain constant and historical depth may be extended mi + 1 = mi + 1, to the benefit of a reduction of (better precision). On the other hand, if Li > 0, true rate ρ is considered to change and the historical depth needs to be reduced mi + 1 = mi − Li, to the benefit of (better accuracy).
At every elementary time step ∆t, the retained count rate estimate is therefore calculated over an adaptable temporal depth mi, Eq. (17) becoming:
Remaining high-frequency fluctuations can now be reduced using a second, recursive smoother, for instance a Brown’s double exponential filter . A first exponential smoothing is performed on with a smoothing parameter δi such as:
A last exponential smoothing is eventually performed on under the form:
The parameter δi changes as a function of the parameter mi and its strength is set with the parameter W:
Finally, the Brown’s estimation is calculated such as:
Figures 5 and 6 illustrate the advantage of the hereby described nonlinear filters over conventional, moving average filters with a 20% rate variation, respectively, in a low count rate configuration (5 counts per sample) and in a higher count rate configuration (500 counts per sample). The nonlinear filter has been set with parameters Qγ = 1.645 (γ = 90%) and W = 0.2 and compared to moving average filters set with m = 50 samples (soft) and m = 500 samples (hard). Figure 5 shows that nonlinear filtering offers a better compromise between precision and accuracy, though the detection of small changes within large statistical fluctuation remains unreachable . At higher count rates (Figure 6), nonlinear filtering permits the detection of the rate change and ensures a significant gain, operating both faster and more precisely than both moving averages.
Such nonlinear smoothing algorithms, easily embedded into programmable components, have for instance been implemented into a Geiger-Müller dosimeter fixed on a wireless robot used for radiological threat detection . This algorithmic building block is plays a key role in the nuclear counting methods studied in the next sections, namely compensation measurements and sensor network processing.
4. Compensation measurement
In many cases, radiation monitoring requires the counting of a signal from a first radiation source within an interference signal induced by a second particle emitter, namely α/β vs. γ; n vs. γ, γ vs. γ … The most efficient techniques consist in the recognition of the particle origin associated with each individual pulse event by coincidence/anti-coincidence, pulse height discrimination (PHD) or pulse shape discrimination (PSD) . However, event-by-event discrimination techniques may be found unreliable in particular mixed field configurations. Compensation methods are an alternative solution when addressing such limitations [17, 18]: their principle lies within measuring count rates ρA from a first detector A, sensitive to all particles, and comparing the result with count rates ρB from a second detector B, only significantly sensitive to background contributions (typically gamma rays). The estimation ρC of the count rate associated with particles of interest is obtained by subtraction of ρA with ρB such as:
where ω is a correction factor taking in account the fact that detector B is not strictly equivalent to detector A in terms of response as a function of the energy and spatial localization of incident background particles. Three challenges are to be faced in compensation measurement:
increase in fluctuation level;
apparition of negative count rates without physical sense;
loss of reliability (impact of energy and anisotropy of the background signal).
Values of ρA, i∆t and ρB, i∆t at the time ti are described by Poisson processes, as already stated in the third section of this chapter. Therefore, if ω = 1, values of ρC, i∆t are described by a Skellam process such as:
The expectation and the variance of the random variable are, respectively, and under the same assumption as in Section 4. The variance definition highlights an increase of fluctuation level in comparison with single-channel measurement. It is therefore required to reduce this variance using a suitable smoothing filter, such as the CST a nonlinear filter described in the previous section (cf. Eqs. (36) and (37)):
Reduced variances and calculated according to Eq. (36) are used to determine as:
If the compensation factor ω remains constant but different from 1, Eq. (42) becomes:
In practice, the factor ω is not constant, due to the impact of energy and spatial distributions of incident particles. We then introduce notations and σ2(ω) for the expectation and variance of the variable ω. A resulting variance is therefore calculated by taking into account both statistical error and bias such as:
The estimation of and σ2(ω) is complicated by the experimental dependence of these parameters. An approach is proposed in [19, 20], in which a database is built from measures acquired in representative areas and in absence of the signal particle of interest: ρC = 0. In these conditions, compensation factors are obtained for each measurement point (1 ≤ q ≤ nq), allowing for an empirical mean and variance to be estimated.
Based on the generalized variance expressed in Eq. (43), a hypothesis test is built to select positive and significant values of . Algorithm 1 presents the detection test in which the presence of particles of interest is detected in compliance with a confidence level γ governing the test. Most of the time, ωq values can be considered to follow a Normal law, which allows us to apply an envelope coverage factor Qγ associated with a confident level γ as:
Then (detection hypothesis H1 is accepted)
Else (detection hypothesis H1 is rejected)
Figure 7 synthetizes the principle, inputs, and outputs of the compensation technique.
The method improves the reliability of compensation measurement with the use of a recorded database. Moreover, accuracy and decision threshold associated with the particles of interest are optimized using an adaptive filter, smoothing individual channels while suppressing all negative or non-significant values. The approach described in the present section has been successfully implemented in varied applications, such as α/β contamination meters or gadolinium-based neutron detectors [20, 21].
As a perspective, it has been demonstrated that the multiplication of channels, such as illustrated in Figure 8, allows the system to learn a prior distribution for the signal over a set of pixels as a function of incident energy and spatial origin of background particles. Dispatching and data along (X > 1) dimensions induces a reduction of detection threshold and thus an improvement of measurement reliability. Compensation factors and variances σ2(ωj, k) are determined for every 1 ≤ j ≤ X A-type detector, and every 1 ≤ k ≤ X B-type detector. Resulting count rates ρC, i are estimated as:
5. Moving source detection
Radiation portal monitors (RPM) are implemented to detect radioactive sources, carried by a vehicle in motion, through the monitoring of a count rate measured by large-volume detectors. Two main issues arise in RPM development: correcting the shadow shielding effect observed when the vehicle is dense enough to impact the baseline of the signal, and improving the detection capability (increasing true detection minus false alarm detection probability).
RPM detection strategy is based on a hypothesis test where the estimated signal at the time ti is continuously compared to a threshold h, itself determined in comparison with the signal distribution under H0. Let θ0 be the expected background count rate without any vehicle in the environment surrounding the RPM. A decision threshold (DT) is set, following the same philosophy as presented previously (Eqs. (20)–(26)), as a function of variance σ2(θ0) and a confidence level associated with a false detection risk α
During the passage of a dense vehicle, θ0 will decrease due to gamma-ray attenuation as the vehicle acts as a radiation shield. Such baseline alteration is noted ωθ0, where ω ϵ [0; 1] is the attenuation factor. An added count rate from a source with intrinsic rate θ1, put onboard the vehicle, will thus lead to a total signal θT = ωθ0 + θ1. If ω = 1 and θ1 > DT, the source is detected with a non-detection risk:
If ω < 1, Eq. (48) becomes:
This effect, so-called “shadow effect,” induces a significant loss in detection capability (β↗) even for ω ≈ 1.
Many works have been done in order to restore the baseline (ω → 1) and all of them use a database recorded when a representative sample of void vehicle in passing through the RPM [22, 23]. An alternative method based on time series analysis has been developed to restore the baseline without using any prior knowledge about the vehicle and the experimental conditions hoping for gain in flexibility . The latter is described below.
In the first place, the minimization of DT requires the implementation of an efficient smoothing filter, minimizing the high-frequency variance σ2(θ0) and subsequently the β risk, while preserving the temporal shape of the signal of interest θ1(ti). Thus, the nonlinear filter CST (Eqs. (36) and (37)) has proven efficient for this purpose. The single-channel RPM estimates the random variable ρi∆t∼(θT∆t) at each time ti such as:
Estimations are continuously recorded into an historical memory with a depth m allowing calculating, at each time step ti, the filtered logarithmic derivative of the signal:
with α1 and l being, respectively, a smoothing parameter and the derivative depth.
The trend of the signal, which can be constant, decreasing or increasing, is represented by a slope state Di with values between −1 and 1 such as:
with α2 > 0 being a parameter for variation significant.
The state of the signal Si is labeled by a number between 1 and 8, defined as illustrated in Figure 9. The first line describes the passage of a vehicle containing a source without shadow effect; the second line corresponds to the passage of a dense vehicle with no source; and the third one to the passage of a dense vehicle containing a radioactive source (shadow shielding).
States Si can be determined with knowledge of Di and Si − 1 using a sequential logic algorithm detailed under the form of a state diagram in Figure 10. To solve the problem, states 3 and 8 automatically pass to state 1 after a preset watchdog time τw.
Knowing the state of the system, the baseline of the signal can be restored. The upper level (UL) (state 1) and the lower level LL (state 6) are firstly estimated in a recorded time series at time ti:
The baseline is restored to obtain corrected count rate estimations such as:
∀k ∈ ⟦1; ξ⟧,
Figure 11 illustrates the baseline restoration: the correction algorithm enables the detection of a source originally hidden by shadow effect. A simulation study has shown the significant gain in detection probability with the maintaining of a stable false detection rate .
The conception of a RPM primarily consists in designing detection blocks with a maximized sensitivity according to the application view and cost-effectiveness strategies. Signal processing is then to be implemented in the system in order to tune its detection capabilities. The improvement of RPM performance forms an active topic of research. It has notably been established that the spectral analysis of the signal, even for unresolved detectors, allows a gain in detection performance . Another upgrade can be achieved by time series analysis techniques, especially when RPM are deployed in a network, which allows the implementation of correlation methods . Figure 12 presents the schematic of a RPM network implementing n channels and dedicated to moving source detection.
The network configuration enables two complementary types of detection: the first one based on traditional temporal analysis of individual channel , the second one based on frequency analysis, searching for a phase ϕ maximizing the correlation between channels . When the network is linear and the source carrier has a constant velocity, this phase corresponds to the periodic echo of the signal increase on the first channel as seen on the other channels. The difference in nature between both methods introduces a quantitative information gain, and thus a potential improvement in detection capability [27–29].
A correlation vector Rφ is calculated by scanning the product of all channels with phase φ such as:
For the vehicle passing from detectors 1 to n.
The algorithm firstly determines a phase φ0 maximizing Rφ, then the significance of the associated temporal correlation is evaluated with a hypothesis test in which H0 is the null hypothesis (no echo detected) and H1 is the detection hypothesis (echo detected). Values of Rφ are compared to the mean and variance of their distribution. Mean and empirical variance σ2(Rφ) are calculated according to:
The detection test reads:
Then H1 is accepted
Else H0 is accepted
The use of the empirical variance σ2(Rϕ) ensures a significant gain in detection capability under challenging signal-to-noise ratios [27, 28]. However, an algorithmic refinement is achieved with the introduction of a Normal law a priori on count rate distributions . Thus, a modified variance σ2(Rϕ) is obtained by the estimation of individual variances σ2(ρj, i) provides by the nonlinear filter (cf. Eq. (37)) for every detector j and memory slot i. Its calculation is presented in the following recursive formula:
The detection algorithm mixes the detection according to each individual channel with the detection using the correlation factor. In both cases, a decision threshold (DT) is calculated as a function of a false detection risk α. Let DTj, ϕ the decision threshold associated with the channel j and the phase ϕ reads:
And let DTRϕ be the decision threshold associated with the correlation factor Rϕ such as
Algorithm 2 presents the mechanism of the cumulative detection.
And if DTRφ − Rφ ≤ 0,
Then, the detection hypothesis H0 is accepted and the hypothesis H1 is rejected,
Else if ,
Then, the detection hypothesis H1 is accepted and the hypothesis H0 is rejected,
Or if ,
Then, the detection hypothesis H1 is accepted and the hypothesis H0 is rejected,
And, the velocity of the source is equal to where L is the distance between detectors.
It has been proven in [27–29], the largely significant added-value in term of detection capability permits by the implementation of the correlation based detection. The true detection rate is increased while maintaining very low false alarm rate. Figure 13 presents a system realized by the CEA which implements the correlation method .
Different count rate processing methods have been presented in this chapter: an adaptive smoother, a background discrimination method and two algorithms improving the detection of moving sources. In these algorithms, frequentist inferences are realized on the basis of measured data. These types of approaches are well suited for real-time processing, allowing taking decision with very few iterations, compared to Bayesian inferences which are more suited for post-processing analyses.
The nonlinear smoother is proved to be a key building block in radiation monitors, delivering a fine estimation of count rate expectation with a minimized associated variance. Both expectation and variance estimations are used to apply hypothesis tests addressing many problematics in radiation monitoring such as for instance those already developed hereby: compensation and RPM network.
- Current ADCs are available with tradeoffs between resolution and sampling frequency such as: 16 bit / 100MS/s and 8 bit / 1 GS/s. CAEN Electronic instrumentation, 724 Digitizer Family, CEAN data sheet, 2015.