## Abstract

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.

### Keywords

- data processing
- nuclear counting
- radiation monitor
- signal processing
- filtering
- frequentist inference

## 1. Introduction

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:

Voltage supply,

Detector part,

Front-end electronics,

User interface.

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 *p*_{X → 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 *p*_{X → Y}(*t*) is usually very small (1/*λ* ≫ *t*). In these conditions, the Binomial law converges toward a Poisson law

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 *t*_{0} = 0 & ∀ *i* ≥ 1, *ti* = *t*_{i − 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: *N*_{i − 1}, *N*_{i − 2}, … recorded in a memory according to the assumption that *p*(*ρi*|*Ni*, *N*_{i − 1}, *N*_{i − 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 [6]:

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

The relative stochastic uncertainty *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 *ti* over the temporal depth *m* + 1 are identically weighted, the estimate *ρ* 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

First approaches consist in the implementation of preset count filter providing a fixed variance *ϑ* is used to assign more weight to recent than older count values such as Eq. (12) become:

Among FIR filters, the exponential moving average (EMA) remains widespread [7, 8], but do not fully deals with the tread-off issue between accuracy and precision.

The algorithm translation of the actualization of *mi* is the building of infinite impulse response (IIR) dedicated to nuclear counting [9]. Such nonlinear filtering requires a hypothesis test to detect the changes in count rate *ρ*. The null hypothesis **H**_{0} and the detection hypothesis **H**_{1} are formalized as follows:

In a first approach [10], 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) [11] and the centered significance test (CST) [12]. 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

In the CST test, the vector *ρ*. For every temporal depth *k* ∈ ⟦1; *mi*⟧, the difference between count rate estimations

The method is based on a comparison between actual and expected distributions of _{0} and H_{1}, respectively [13]. The distribution *℘*

The expectation *i* − *mi*)∆*t* and *i*∆*t*. Moreover, we will make use of assumption

We will note

Under **H**_{0}, *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

In practice, for embedded implementations, it is impossible to sample and interpolate distributions *i* and *k*. Moreover, when *i* and *k*, *Q*_{1 − α} is a quantile of

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 *Q*_{1 − α} such as:

If **H**_{0} is accepted with a confidence level equal to 1 − *α*.

Under **H**_{1}, *DL* is defined as the detection limit. *DL* is determined in compliance with a given risk of non-detection *β* = *p*(**H**_{0}| **H**_{1}). *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 *Q*_{1 − β} of the error function erf such as:

An equivalent confidence level 1 − *α* = 1 − *β* = *γ* is considered, and Eq. (23) is solved recursively such as:

With

When *y* → ∞,

If **H**_{1} is accepted with a confidence level 1 − *β* = *γ*.

The number *Li* of significant changes recorded into memory *m*_{i + 1}:

If *Li* = 0, true rate *ρ* is considered to remain constant and historical depth may be extended *m*_{i + 1} = *mi* + 1, to the benefit of a reduction of *Li* > 0, true rate *ρ* is considered to change and the historical depth needs to be reduced *m*_{i + 1} = *mi* − *Li*, to the benefit of

At every elementary time step ∆*t*, the retained count rate estimate *mi*, Eq. (17) becoming:

With

This nonlinear approach performs advantageously in comparison with conventional linear filters [12, 14], allowing to maintain sufficient precision while rate changes in the signal occur.

Remaining high-frequency fluctuations can now be reduced using a second, recursive smoother, for instance a Brown’s double exponential filter [14]. A first exponential smoothing *δi* such as:

With

A last exponential smoothing

With

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

With

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

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 [15]. 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) [16]. 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

The expectation and the variance of the random variable *CST* a nonlinear filter described in the previous section (cf. Eqs. (36) and (37)):

Reduced variances

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 *σ*^{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 *σ*^{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 *q* ≤ *nq*), allowing for an empirical mean

Based on the generalized variance expressed in Eq. (43), a hypothesis test is built to select positive and significant values of *ω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:

Algorithm 1:

If

Then **H**_{1} is accepted)

Else **H**_{1} 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

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 *X* > 1) dimensions induces a reduction of detection threshold *σ*^{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 *ti* is continuously compared to a threshold *h*, itself determined in comparison with the signal distribution under **H**_{0}. 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 [24]. 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

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 *S*_{i − 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

∀*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 [24].

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 [25]. 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 [26]. 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 *ϕ* maximizing the correlation between channels

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 **H**_{0} is the null hypothesis (no echo detected) and **H**_{1} is the detection hypothesis (echo detected). Values of *Rφ* are compared to the mean and variance of their distribution. Mean *σ*^{2}(*Rφ*) are calculated according to:

The detection test reads:

Algorithm 2:

If

Then **H**_{1} is accepted

Else **H**_{0} 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 [29]. 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:

with

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 *DT*_{j, ϕ} 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.

Algorithm 3:

If

And if *DTRφ* − *Rφ* ≤ 0,

Then, the detection hypothesis H_{0} is accepted and the hypothesis H_{1} is rejected,

Else if

Then, the detection hypothesis H_{1} is accepted and the hypothesis H_{0} is rejected,

Or if

Then, the detection hypothesis H_{1} is accepted and the hypothesis H_{0} is rejected,

And, the velocity of the source is equal to *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 [30].

All of these algorithms will be implemented in a dedicated DSP card [3] and the compliance of the RPM system will the standard ANSI42-35 will be tested in due course [31].

## 6. Conclusion

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.

## Notes

- 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.