Recent Developments in Count Rate Processing Associated with Radiation Monitoring Systems Recent Developments in Count Rate Processing Associated with Radiation Monitoring Systems

This chapter presents some recent data processing developments associated with radia- tion 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.


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 fieldprogrammable 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 N i 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][2][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.

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 0 ) of the decay instant, where t 0 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 = 0 being the start of the observation) is obtained as: In a radioactive source containing a population of N X unstable nuclei, the decay of an individual nucleus does not impact the decay of the others. The Bernoulli trial is therefore repeated N X 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, N X 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 P such as: Expectation and variance of the number of decays are equal to N X λ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 ρ = εηN X λ thus becomes the parameter of the distribution of measured count values before t: At each time t i , sampled such as t 0 = 0 & ∀ i ≥ 1, t i = t i À 1 + Δt, the raw estimation of the count rate ρ(t i ) = ρ i is provided by measuring N i , 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 t i maximizing both precision σ(ρ i ) ! 0 and accuracy σ(t i ) ! 0. Algorithmic techniques to meet this expectation are discussed in the next section.

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 N i À 1 , N i À 2 , … recorded in a memory according to the assumption that p(ρ i |N i , N i À 1 , N i À 2 , …) is more precise than p(ρ i |N i ). 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 σ 2 b ρ i À Á can be estimated as: The relative stochastic uncertainty σ b ρ i À Á =b ρ i 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 b t i effectively corresponding to the current count rate estimate b ρ i . Because the sampling times t i over the temporal depth m + 1 are identically can only be minimized to the detriment of σ b t i , leading to a degradation of accuracy when ρ is varying. One way to address this issue is to actualize the temporal depth m i after every count rate estimation. The optimal value for m i is a function of the temporal behavior of the rate ρ at the time t i : 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: 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 m i 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 b ρ k i to designate both the underlying random variable and its actual values.
In the CST test, the vector b ρ k, 1 ≤ k ≤ mi i is scanned to find a potential change in the true rate ρ. For every temporal depth k ∈ ⟦1; m i ⟧, the difference between count rate estimations Δb is the quantity which will be tested for significance.
The method is based on a comparison between actual and expected distributions of Δb ρ k i under H 0 and H 1 , respectively [13]. The distribution D of Δb ρ k i is the difference between two weighted The expectation E Δb ρ as common in nuclear counting experiments will finite statistics [1]. The variances associated with both uncorrelated random processes are summed to obtain a cumulative standard deviation for Δb ρ k i . According to the equality between expectation and variance, we obtain: We will note D E Δb ρ , the distribution of the difference random variable with its first and second order moments.
is the quantile of the error function (err) with a confidence level 1 À α k i : In practice, for embedded implementations, it is impossible to sample and interpolate distri- Under this assumption, for every value of i and k, quantile of N 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 Q 1 À α such as: where 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 Δb ρ k i ≥ DL k i , hypothesis H 1 is accepted with a confidence level 1 À β = γ.
The number L i of significant changes recorded into memory Δb ρ k i is used to calculate the next value of temporal depth m i + 1 : If L i = 0, true rate ρ is considered to remain constant and historical depth may be extended m i + 1 = m i + 1, to the benefit of a reduction of σ b ρ À Á (better precision). On the other hand, if L i > 0, true rate ρ is considered to change and the historical depth needs to be reduced m i + 1 = m i À L i , to the benefit of σ b t (better accuracy).
At every elementary time step Δt, the retained count rate estimate b ρ * i is therefore calculated over an adaptable temporal depth m i , 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 b ρ 1 i is performed on b ρ * i with a smoothing parameter δ i such as: With The parameter δ i changes as a function of the parameter m i and its strength is set with the parameter W: Finally, the Brown's estimation b ρ * * i is calculated such as: 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  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.

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: Figure 6. Behavior of smoothing filters over a 20% rate variation at 500 counts per sample.
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 t i 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 Sk ρ A, i Δt; ρ B, i Δt such as: The expectation and the variance of the random variable b ρ C, i Δt are, respectively, Reduced variances σ 2 b ρ * * A, i and σ 2 b ρ * * B, i calculated according to Eq. (36) are used to determine σ 2 b ρ * * C, i 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, are obtained for each measurement point (1 ≤ q ≤ n q ), allowing for an empirical mean ω ¼ 1 Based on the generalized variance expressed in Eq. (43), a hypothesis test is built to select positive and significant values of b ρ C, i . 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: 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 Q γ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 2σ 2 b ρ * * C, i r 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 b ρ * * * A, q and b ρ * * * B, q data along (X > 1) dimensions induces a reduction of detection threshold Q γ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 2σ 2 b ρ * * C, i r and thus an improvement of measurement reliability. Compensation factors ω j, k 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:

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 b ρ i at the time t i 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 ω e [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 (t i ). 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$P(θ T Δt) at each time t i such as: Estimations are continuously recorded into an historical memory with a depth m allowing calculating, at each time step t i , the filtered logarithmic derivative_ ρ i 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 D i with values between À1 and 1 such as: with α 2 > 0 being a parameter for variation significant.
The state of the signal S i 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 S i can be determined with knowledge of D i 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 t i : The baseline is restored to obtain corrected count rate estimations b ρ * * * 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 [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 H 1 |b ρ 1 , …, b ρ n , the second one based on frequency analysis, searching for a phase ϕ maximizing the correlation between channels H 1 |φρ 1 , …,ρ n . 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][28][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 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 R and empirical variance σ 2 (R φ ) are calculated according to: Figure 12. Schematic of a system based on correlation detection.
The detection test reads: Algorithm 2: 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: 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 DT R ϕ be the decision threshold associated with the correlation factor R ϕ such as ∀φ ∈ 1; ξ À 1 n À 1 " # ! ! , Algorithm 2 presents the mechanism of the cumulative detection.
Algorithm 3: Then, the detection hypothesis H 0 is accepted and the hypothesis H 1 is rejected, Then, the detection hypothesis H 1 is accepted and the hypothesis H 0 is rejected, 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 φ 0 Δt where L is the distance between detectors.
It has been proven in [27][28][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].

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.