Approximate sequential Bayesian filtering to estimate ²²²Rn emanation from ²²⁶Ra sources using spectral time series

A new approach to assess the emanation of ²²²Rn from ²²⁶Ra sources based on γ-ray spectrometric measurements is presented. While previous methods have resorted to steady-state treatment of the system, the method presented incorporates well-known radioactive decay kinetics into the inference procedure through the formulation of a theoretically motivated system model. The validity of the ²²²Rn emanation estimate is thereby extended to regimes of changing source behavior, potentially enabling the development of source surveillance systems in the future. The inference algorithms are based on approximate recursive Bayesian estimation in a switching linear dynamical system, allowing regimes of changing emanation to be identified from the spectral time series while providing reasonable filtering and smoothing performance in steady-state regimes. The derived method is applied to an empirical γ-ray spectrometric time series obtained over 85 d and is able to provide a time series of emanation estimates consistent with the physics of the emanation process.

1 Introduction

²²²Rn is an odorless, colorless, radioactive noble gas, occurring naturally in the environment as part of the primordial ²³⁸U decay chain. Due to its high mobility in the environment, radon can accumulate inside buildings where, in conjunction with its decay products (short-lived progeny; SLP), it is responsible for the most significant natural exposure of the general public to ionizing radiation, which is an important consideration for incidents of lung cancer (Darby et al., 2005). In the 2013/59/EURATOM treaty, 300 Bq m⁻³ was stipulated as the action level for indoor ²²²Rn activity concentrations, beyond which mitigation is required. On the other hand, however, measurements of outdoor ²²²Rn activity concentration at environmental levels have numerous beneficial applications in environmental sciences, including (but not limited to): as a tracer of terrestrial influence on air masses (Chambers et al., 2016, 2018); as a tool for classifying the atmospheric mixing state (e.g., Perrino et al., 2001; Williams et al., 2013, 2016; Chambers et al., 2019b, a); and as a tool for estimating integrated local- to regional-scale emissions of trace gases with similarly distributed sources, such as CH₄, N₂O, or CO₂ (Levin, 1987; Biraud et al., 2002; Laan et al., 2010; Levin et al., 2021). For these reasons, there is an interdisciplinary need for SI-traceable calibration procedures at low activity concentrations for atmospheric radon monitors, and the associated realization and dissemination of the unit Bq m⁻³, which can only be feasibly realized through emanation sources rather than gaseous standards (Mertes et al., 2020; Linzmaier and Röttger, 2013; Röttger et al., 2014). In these prior studies, a method was presented that enabled γ-ray spectrometric data from an open ²²⁶Ra source (i.e. emitting ²²²Rn) to be used to estimate the resulting activity concentration of ²²²Rn in a closed volume, by measuring the activity of residual SLP in the source to quantify the ²²²Rn emanation. In recent years, there has been a trend towards the use of dynamic calibration conditions for ²²²Rn (e.g., Fialova et al., 2020), eliminating the need for the long build-up period associated with static conditions. However, it has been demonstrated that emanation from most materials demonstrates a strong dependence on humidity and temperature as a result of changes in diffusion properties (Janik et al., 2015; Stranden et al., 1984; Zhou et al., 2020), which generally results in a correlation between emanation and environmental parameters. In these cases, previously established methods to determine the amount of emanating ²²²Rn fail over a considerable time span, since the dynamical processes taking place in the source are not accounted for. Hence, experimental investigations of source behavior under different environmental conditions are hardly possible using established methods. Yet, this capability is particularly important for the realization of in situ field calibrations of large volume atmospheric ²²²Rn monitors, since dynamic methods would simplify the technical aspects of in situ field calibrations considerably. Said limitations are pointed out and discussed in the theoretical section of this work and have not been stated nor addressed elsewhere. A possible way to correct for such environmental influences, however, lies in determining the amount of emanating ²²²Rn in near real time. We present herein that this can be achieved through continuous measurement of spectrometric time series of the ²²²Rn emanation sources and a suitable method for data analysis that addresses the dynamic behavior of the system. The main contributions of the present work are the discussion of the limitations of established methods and the derivation and implementation of an alternative method, based on the well established computational techniques of recursive Bayesian estimation. First results obtained by application of the proposed method to experimental data are presented, which illustrate the theoretically motivated limitations and which can be well explained on the basis of the physical processes taking place in the emanation source, justifying the correctness and superiority of the presented method.

2 Theory and derivations

2.1 Radioactive decay kinetics

Radioactive decay chain kinetics are described by a linear time invariant (LTI) system of ordinary differential equations. Historically, this has been expressed in terms of the Bateman equations (Bateman, 1910), but recently this has been more conveniently written in matrix form (Pressyanov, 2002; Levy, 2019; Amaku et al., 2010). Undistorted radioactive decay kinetics of a decay chain of n nuclides without, or with only negligible, branching, as in the case of the ²²⁶Ra decay chain, can thus be written concisely as Eq. (1). The fundamental matrix K can be constructed such that it consists of the respective decay constants λ_i on its diagonal and generally on its first superdiagonal, while all other entries are 0. A denotes a vector consisting of the activities of the respective nuclides in the decay chain. Equation (1) can readily be discretized using the matrix exponential of K, which is in a form that can be conveniently computed by diagonalization of K, e.g., through its (symbolically accessible) Jordan canonical form or its Eigen decomposition.

$$\begin{array}{}\text{(1)}& \begin{array}{rl}\mathrm{d}\mathbf{A}& =\mathrm{d}\left[\begin{array}{c}{A}_{\mathrm{0}}\\ A_{\mathrm{1}}\\ \mathrm{\dots }\\ A_n\end{array}\right]\\ & =\left[\begin{array}{cccc}-{\mathit{\lambda }}_{\mathrm{0}}& {\mathit{\lambda }}_{\mathrm{0}}& \mathrm{0}& \mathrm{⋮}\\ \mathrm{0}& -{\mathit{\lambda }}_{\mathrm{1}}& {\mathit{\lambda }}_{\mathrm{1}}& \mathrm{⋮}\\ \mathrm{\cdots }& \mathrm{\cdots }& \mathrm{\cdots }& \mathrm{⋮}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& -{\mathit{\lambda }}_n\end{array}\right]\left[\begin{array}{c}{A}_{\mathrm{0}}\\ A_{\mathrm{1}}\\ \mathrm{\dots }\\ A_n\end{array}\right]\mathrm{d}t\\ & =\mathbf{KA}\mathrm{d}t\end{array}\end{array}$$

2.2 ²²²Rn emanation

The release of ²²²Rn from a ²²⁶Ra source distorts the dynamics described above, due to the introduction of an additional sink-term η (²²²Rn atoms released per unit time), which is not directly quantifiable experimentally. Previously, a method was presented to measure a steady-state emanation coefficient, previously understood to be the ratio of exhaled and generated ²²²Rn atoms (Mertes et al., 2020; Linzmaier and Röttger, 2013) at any instant in time.

$$\begin{array}{}\text{(2)}& \mathit{\chi }={\displaystyle \frac{\mathit{\eta }}{A_{\text{Ra-226}}^{\mathrm{S}}}},\end{array}$$

where χ is the emanation coefficient, η is the release of ²²²Rn atoms per unit time, and $A_{\text{Ra-226}}^{\mathrm{S}}$ is the ²²⁶Ra activity of the source.

By first principles, however, the ²²²Rn activity $A_{\text{Rn-222}}^{\mathrm{S}}$ must follow first-order continuity as in Eq. (3).

$$\begin{array}{}\text{(3)}& \mathrm{d}{A}_{\text{Rn-222}}^{\mathrm{S}}=-{\mathit{\lambda }}_{\text{Rn-222}}{A}_{\text{Rn-222}}^{\mathrm{S}}\mathrm{d}t+{\mathit{\lambda }}_{\text{Rn-222}}{A}_{\mathrm{Ra}-\mathrm{226}}^{\mathrm{S}}\mathrm{d}t-{\mathit{\lambda }}_{\text{Rn-222}}\mathit{\eta }\mathrm{d}t,\end{array}$$

where $A_{\text{Rn-222}}^{\mathrm{S}}$ and $A_{\mathrm{Ra}-\mathrm{226}}^{\mathrm{S}}$ are the activities of the respective nuclides in the source material.

To this point, the best method to derive χ has been by measuring the ratio of SLP and ²²⁶Ra activities within the source, most commonly by γ-ray spectrometry, where χ is defined as

$$\begin{array}{}\text{(4)}& \mathit{\chi }=\mathrm{1}-{\displaystyle \frac{A_{\text{Rn-222}}^{\mathrm{S}}}{A_{\mathrm{Ra}-\mathrm{226}}^{\mathrm{S}}}}=\mathrm{1}-{\displaystyle \frac{A_{\mathrm{Pb}-\mathrm{214}}^{\mathrm{S}}}{A_{\mathrm{Ra}-\mathrm{226}}^{\mathrm{S}}}}.\end{array}$$

According to Eqs. (2) and (3) however, this is only applicable under steady-state conditions. This limitation was not acknowledged or discussed in earlier work. Consequently, application of this simple method is restricted either to regimes of a completely stabilized source, or a closed accumulation volume for the emitted ²²²Rn. In closed volumes, the errors associated with neglecting the dynamics in the source are negligible, since the volumetric ²²²Rn follows the same dynamics. As such, after initial equilibration, Eq. (4) holds, even with changing χ. However, in this case χ is more accurately understood to be a partitioning coefficient of ²²²Rn between the source and the closed volume, rather than the definition given in Eq. (2). The result of Eq. (3) is that the measurable $A_{\text{Rn-222}}^{\mathrm{S}}\left(t\right)$ is given by the convolution of the latent η(t) with an impulse response function that is defined by the radioactive decay kinetics. Thus, the estimation of η(t), or equivalently χ(t), based on measurements of $A_{\text{Rn-222}}^{\mathrm{S}}\left(t\right)$ is an inverse problem and cannot be carried out feasibly by simple numerical estimation of the gradient (cf. Eq. 3) due to the ubiquitous Poisson noise associated with radioactivity measurements. Moreover, in dynamic calibrations, η(t) must be expected to vary with changes in the environmental conditions. It is also expected that this dependency will be strongly related to the source design and its specific production parameters. A further consideration is that, when using this method, it is not readily possible to accurately measure correlations of χ(t) with environmental conditions on timescales smaller than at least five half-lives of ²²²Rn without considering the decay kinetics. These limitations make it infeasible to use the existing direct approach to continuously estimate ²²²Rn release as required in dynamic calibration procedures, or to derive correction factors for different environmental conditions, since the time required would be far too long for such measurements due to the half-life of ²²²Rn of approximately 3.8 d. Here we present and discuss a new approach that embraces the described behavior and enables the release of ²²²Rn from sources to be estimated more accurately using continuous spectrometric measurements with the generalization to non-steady-state situations. The algorithms and assumptions presented have been chosen such that, in the future, the necessary computations would be feasible on relatively low-power devices (e.g., current single-board computers) in an online fashion.

2.3 Recursive Bayesian estimation and model formulation

Recursive Bayesian estimation describes a class of algorithms to perform statistical inference in dynamical systems that can be modeled by a (first-order) Markov process. The general idea of these methods is to sequentially form priors for a state vector x and a dynamical model of the system, and use noisy measurements y related to x to correct them through a measurement model and Bayes theorem (Särkkä, 2013). This method can be used to perform statistical inversion, in that it enables the estimation of latent terms whose values in a specific dynamical system are only partly measured. This is closely related to the case described in the previous section, given that the state vector, the dynamics, and the measurements can all be suitably modeled. In this setting, two collections of conditional probability distributions are of interest, which are commonly called the filtering distributions $p\left({\mathit{x}}_n|{\mathit{y}}_{\mathrm{1}:n}\right)$ and the smoothing distributions $p\left({\mathit{x}}_n|{\mathit{y}}_{\mathrm{1}:N}\right)$, where the notations 1:n and 1:N denote the collection of all data observed up to time t_n and the collection of all data, respectively, and t specifies an instant in time in the set of measurement times, T, indexed by n. In cases where x follows a first-order Markov process, and under certain conditional independence assumptions, the definition of $p\left({\mathit{x}}_n|{\mathit{y}}_{\mathrm{1}:n}\right)$ and $p\left({\mathit{x}}_n|{\mathit{y}}_{\mathrm{1}:N}\right)$ can be expressed recursively (Särkkä, 2013; Särkkä and Solin, 2019). Prediction of the state vector x_n at time t_n given a collection of measurements ${\mathit{y}}_{\mathrm{1}:n-\mathrm{1}}$ up to time t_n−1 and the filtering density of the state at time t_n−1, $p\left({\mathit{x}}_{n-\mathrm{1}}|{\mathit{y}}_{\mathrm{1}:n-\mathrm{1}}\right)$ is given by the Chapman–Kolmogorov Eq. (5) (Särkkä and Solin, 2019).

$$\begin{array}{}\text{(5)}& p\left({\mathit{x}}_n|{\mathit{y}}_{\mathrm{1}:n-\mathrm{1}}\right)=\int p\left({\mathit{x}}_n|{\mathit{x}}_{n-\mathrm{1}}\right)p\left({\mathit{x}}_{n-\mathrm{1}}|{\mathit{y}}_{\mathrm{1}:n-\mathrm{1}}\right)\mathrm{d}{\mathit{x}}_{n-\mathrm{1}}\end{array}$$

Upon observation of y_n, the density predicted by Eq. (5) is corrected into a filtered posterior using Bayes theorem with the measurement likelihood p(y_n|x_n), assumed to be conditionally independent of the past states and measurements, i.e. $p\left({\mathit{y}}_n|{\mathit{x}}_{\mathrm{1}:n},{\mathit{y}}_{\mathrm{1}:n-\mathrm{1}}\right)=p\left({\mathit{y}}_n|{\mathit{x}}_n\right)$ (Särkkä, 2013).

Conversely, smoothing refers to computing the density of the state given all measurements in a specific time interval, or the complete collection. In most cases, smoothing can be defined recursively using information inferred from the filtering and starting a backward recursion at the last time instant at which the smoothing and filtering densities are equal. Formally, the smoothing density is recursively defined by Eq. (6) for the above conditional independence assumptions.

$$\begin{array}{}\text{(6)}& \begin{array}{rl}p& \left({\mathit{x}}_n|{\mathit{y}}_{\mathrm{1}:N}\right)\\ & =p\left({\mathit{x}}_n|{\mathit{y}}_{\mathrm{1}:n}\right)\int {\displaystyle \frac{p\left({\mathit{x}}_{n+\mathrm{1}}|{\mathit{x}}_n\right)}{p\left({\mathit{x}}_{n+\mathrm{1}}|{\mathit{y}}_{\mathrm{1}:n}\right)}}p\left({\mathit{x}}_{n+\mathrm{1}}|{\mathit{y}}_{\mathrm{1}:N}\right)\mathrm{d}{\mathit{x}}_{n+\mathrm{1}}\end{array}\end{array}$$

The most notable examples of these types of algorithms are the Kalman filter (Kalman, 1960) and the Rauch–Tung–Striebel smoother (Rauch et al., 1965), which yield the optimal estimators for discrete systems of linear dynamics and independent Gaussian noise, where the filtering and smoothing can be carried out in O(T) time.

For the problem at hand, a stochastic differential equation (SDE) is needed to express the (time-varying) uncertainty associated with the latent continuous variable η(t). This can be seen as an application of the latent-force models introduced in Alvarez et al. (2009), whose link with Bayesian filtering was previously established in Särkkä et al. (2019) and Hartikainen and Särkkä (2012). The specific choice of SDE is subjective but allows the properties of the resultant functions to be constrained, and can be used to encode prior knowledge (Särkkä and Solin, 2019). In this work, without any claims of optimality, the choice was made to use the zero-mean, mean-reverting Ornstein–Uhlenbeck process for the first derivative of η, i.e., Eq. (7), resulting in constraining η to somewhat smooth functions of certain autocorrelation.

$$\begin{array}{}\text{(7)}& {\displaystyle \frac{{\mathrm{d}}^{\mathrm{2}}\mathit{\eta }}{\mathrm{d}{t}^{\mathrm{2}}}}=-\mathit{\gamma }{\displaystyle \frac{\mathrm{d}\mathit{\eta }}{\mathrm{d}t}}+\mathit{\sigma }\mathrm{d}{\mathit{\beta }}_t,\end{array}$$

where dβ_t describes the increments of a scalar standard Brownian motion.

While this choice is not entirely representative of the physical mechanisms related to η, in that, for example, it allows for negative values of η, it results in a Gaussian process for which the inference procedure has convenient conjugacy and thus a closed-form solution, such that the resultant algorithms are suitable for online operation on low-power computational hardware that can be reasonably used to monitor an emanation source during its operation. It should also be noted that Eq. (2) was formulated with this in mind, in the sense that η is modeled as state-independent, rather than state-dependent, as would likely be more realistic for a diffusive process like ²²²Rn emanation. However, these theoretical inaccuracies did not manifest in practice with the experimental data presented in Sect. 3, given that η is far enough from 0 and the collected data are strong enough, while only approximate inference is of interest. The system is thus modeled to follow the combined SDE given in Eqs. (8)–(11) in terms of the Itō stochastic integral (Särkkä and Solin, 2019; Hartikainen and Särkkä, 2012).

$$\begin{array}{}\text{(8)}& {\displaystyle }\mathrm{d}\mathit{x}=\mathbf{K}\mathit{x}\mathrm{d}t+\mathbf{L}\mathrm{d}{\mathit{\beta }}_t\text{(9)}& {\displaystyle }\begin{array}{rl}& \mathit{x}=\left[\begin{array}{c}{A}_{\text{Rn-222}}^{\mathrm{S}}\\ A_{\mathrm{Ra}-\mathrm{226}}^{\mathrm{S}}\\ \mathit{\eta }\\ \frac{\mathrm{d}\mathit{\eta }}{\mathrm{d}t}\end{array}\right];\\ & \mathbf{K}=\left[\begin{array}{cccc}-{\mathit{\lambda }}_{\text{Rn-222}}& {\mathit{\lambda }}_{\text{Rn-222}}& -{\mathit{\lambda }}_{\text{Rn-222}}& \mathrm{0}\\ \mathrm{0}& -{\mathit{\lambda }}_{\mathrm{Ra}-\mathrm{226}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& -\mathit{\gamma }\end{array}\right];\mathbf{L}=\left[\begin{array}{c}\mathrm{0}\\ \mathrm{0}\\ \mathrm{0}\\ \mathit{\sigma }\end{array}\right]\end{array}\text{(10)}& {\displaystyle }{\mathit{x}}_t={\mathrm{e}}^{\mathbf{K}(t-t_{\mathrm{0}})}{\mathit{x}}_{\mathrm{0}}+\underset{t_{\mathrm{0}}}{\overset{t}{\int }}{\mathrm{e}}^{\mathbf{K}(t-\mathit{\tau })}\mathbf{L}\mathrm{d}{\mathit{\beta }}_{\mathit{\tau }}\text{(11)}& {\displaystyle }p\left({\mathit{x}}_n\right|{\mathit{x}}_{n-\mathrm{1}})\propto \mathcal{N}\left({\mathrm{e}}^{\mathbf{K}(t_n-t_{n-\mathrm{1}})}{\mathit{x}}_{n-\mathrm{1}},\right.\underset{t_{n-\mathrm{1}}}{\overset{t_n}{\int }}{\mathrm{e}}^{\mathbf{K}(t_n-\mathit{\tau })}{\mathbf{LL}}^{\mathbf{T}}{\mathrm{e}}^{{\mathbf{K}}^{\mathbf{T}}(t_n-\mathit{\tau })}\mathrm{d}\mathit{\tau })\end{array}$$

2.4 Measurement model for integrating spectrometric data

Unlike most applications of such filtering algorithms for discretized LTI systems, here the supporting measurements cannot be made at instantaneous moments in time because disintegrations of a specific nuclide can only be recorded over a finite time interval, r_n, indexed by n. In routine radioactivity analysis, this behavior results in decay or ingrowth during measurement corrections. However, in the present case, an additional contribution to the uncertainty results from unknown changes of η over the integration time.

The way we have chosen to model this behavior in the present study is to start by stating that the uncorrupted (i.e., noise-free) measurements are given by Eq. (12), where, for now, we assume that H is known deterministically.

$$\begin{array}{}\text{(12)}& \mathit{y}\left(t,r\right)=\mathbf{H}\underset{\mathrm{0}}{\overset{r}{\int }}\mathit{x}\left(t+\mathit{\tau }\right)\mathrm{d}\mathit{\tau },\end{array}$$

where H is a matrix that maps the state integral to the measurement space.

The elements of H are related to the counting efficiency of the setup and/or nuclide. Note that, in principle, it would be possible to choose H to directly model some region of the spectrum, but we chose to use derived peak areas or even spectrum integrals as the input data such that the elements of H are just the counting efficiencies. More elaborate modeling in the spectrum space was not attempted, since the filtering algorithms generally require inversion of the residual covariance matrix in the measurement space, and thus scale approximately with O(k³) where k is the dimensionality of the measurements.

Since integration is a linear operator, the joint distribution of x_n−1, x_n, and y_n clearly has a Gaussian density, and hence $p\left({\mathit{x}}_n,{\mathit{y}}_n|{\mathit{y}}_{\mathrm{1}:n-\mathrm{1}}\right)$ is readily found by marginalizing over x_n−1, given that the previous time step posterior filtering distribution $p\left({\mathit{x}}_{n-\mathrm{1}}|{\mathit{y}}_{\mathrm{1}:n-\mathrm{1}}\right)$ is already known (and Gaussian). This joint density is derived in Appendix A1. It is assumed that the measurement y_n is related to the state x_n by integrating from $t_{n-\mathrm{1}}+{\mathit{\delta }}_n$ to $t_n=t_{n-\mathrm{1}}+{\mathit{\delta }}_n+r_n$, such that arbitrary integration intervals r_n and arbitrary time offsets δ_n are possible (e.g., if measurements are skipped or delayed). In other words, each time instant where the density of x is inferred coincides with the endpoint of each spectrum acquisition. This approach intrinsically accounts for the ingrowth and decay during the finite integration time, but more importantly, additionally estimates the uncertainty arising from the possible change of η within the integration time. Given that $p\left({\mathit{x}}_n,{\mathit{y}}_n|{\mathit{y}}_{\mathrm{1}:n-\mathrm{1}}\right)$ under the present model is thus accessible, the posterior filtering distribution $p\left({\mathit{x}}_n|{\mathit{y}}_{\mathrm{1}:n}\right)$ is found by conditioning onto the observed value for y_n using the well-known conditioning formula for multivariate Gaussians.

Radioactivity measurements generally follow Poisson statistics. In the framework of Bayesian inference and recursive Bayesian estimation, non-Gaussian noise considerably complicates the inference procedure, since the measurements are then no longer a Gaussian process and thus no longer conjugate with the state. For this reason, there is no exact closed-form solution for the filtering recursions in the case of non-Gaussian noise. Considerable work has been done to address this, including sampling procedures like Markov chain Monte Carlo, particle filtering, or expectation propagation (Minka, 2013), or by assuming that all arising probability density functions (PDFs) are Gaussian combined with an estimation strategy for the moments (e.g., unscented Kalman filter; Julier et al., 2000). Generally, we found such approaches unsuitable due to their computational complexity but also because the measurements can be well approximated as Gaussian (due to the high number of counts that are being observed). Instead, we approximate the measurement likelihood as Gaussian where the moments are evaluated from the previous time step filtering distribution of x. Ebeigbe et al. (2020) proposed that the covariance matrix could be estimated from the mean of the filtered state in a different context, which we adopted. Combining this with the results of Appendix A1, we obtain Eq. (13) as an approximation of the joint density $p\left({\mathbf{x}}_n,{\mathbf{y}}_n\mathrm{|}{\mathbf{y}}_{\mathrm{1}:n-\mathrm{1}}\right)$.

$$\begin{array}{}\text{(13)}& p\left({\mathit{x}}_n,{\mathit{y}}_n\mathrm{|}{\mathit{y}}_{\mathrm{1}:n-\mathrm{1}}\right)\propto \mathcal{N}\left(\left[\begin{array}{c}{\mathbf{F}}_r{\mathbf{F}}_{\mathit{\delta }}{\mathit{\mu }}_{n-\mathrm{1}}\\ {\mathbf{HMF}}_{\mathit{\delta }}{\mathit{\mu }}_{n-\mathrm{1}}\end{array}\right],\left[\begin{array}{cc}{\mathbf{\Sigma }}^{\mathit{x}\mathit{x}}& {\mathbf{\Sigma }}^{\mathit{x}\mathit{y}}\\ {\mathbf{\Sigma }}^{\mathit{x}{\mathit{y}}^{\mathbf{T}}}& {\mathbf{\Sigma }}^{\mathit{y}\mathit{y}}\end{array}\right]\right),\end{array}$$

where

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathbf{\Sigma }}^{\mathit{x}\mathit{x}}={\mathbf{F}}_r\left({\mathbf{F}}_{\mathit{\delta }}{\mathbf{\Sigma }}_{n-\mathrm{1}}{\mathbf{F}}_{\mathit{\delta }}^{\mathbf{T}}+{\mathbf{U}}_{\mathit{\delta }}\right){\mathbf{F}}_r^{\mathbf{T}}+{\mathbf{U}}_r,\\ {\displaystyle }& {\displaystyle }{\mathbf{\Sigma }}^{\mathit{x}\mathit{y}}={\mathbf{F}}_r\left({\mathbf{F}}_{\mathit{\delta }}{\mathbf{\Sigma }}_{n-\mathrm{1}}{\mathbf{F}}_{\mathit{\delta }}^{\mathbf{T}}+{\mathbf{U}}_{\mathit{\delta }}\right){\mathbf{M}}^{\mathbf{T}}{\mathbf{H}}^{\mathbf{T}}+{\mathbf{CH}}^{\mathbf{T}},\\ {\displaystyle }& {\displaystyle }{\mathbf{\Sigma }}^{\mathit{y}\mathit{y}}=\mathbf{HM}\left({\mathbf{F}}_{\mathit{\delta }}{\mathbf{\Sigma }}_{n-\mathrm{1}}{\mathbf{F}}_{\mathit{\delta }}^{\mathbf{T}}+{\mathbf{U}}_{\mathit{\delta }}\right){\mathbf{M}}^{\mathbf{T}}{\mathbf{H}}^{\mathbf{T}}+{\mathbf{HBH}}^{\mathbf{T}}+\mathbf{O}+\mathbf{R}\end{array}$$

and where the index n has been dropped on r and δ for notational brevity and the included matrices are given as follows: R is a variance term that accounts for the variance contribution of the background count rate, which is estimated from prior background measurements and is assumed constant over time, and μ_n−1 and Σ_n−1 denote the mean and covariance matrix of the filtering distribution at the previous time step.

$$\begin{array}{ll}{\displaystyle }{\mathbf{F}}_a=& {\displaystyle }{e}^{\mathbf{K}a}{\mathbf{U}}_a=\underset{\mathrm{0}}{\overset{a}{\int }}{e}^{\mathbf{K}\left(a-\mathit{\tau }\right)}{\mathbf{LL}}^{\mathbf{T}}{e}^{{\mathbf{K}}^{\mathbf{T}}\left(a-\mathit{\tau }\right)}\mathrm{d}\mathit{\tau }\\ {\displaystyle }{\mathbf{M}}_r=& {\displaystyle }\underset{\mathrm{0}}{\overset{r}{\int }}{e}^{\mathbf{K}\mathit{\tau }}\mathrm{d}\mathit{\tau }{\mathbf{C}}_r=\underset{\mathrm{0}}{\overset{r}{\int }}\underset{\mathit{\tau }}{\overset{r}{\int }}{e}^{\mathbf{K}\left(r-\mathit{\tau }\right)}{\mathbf{LL}}^{\mathbf{T}}{e}^{{\mathbf{K}}^{\mathbf{T}}\left(a-\mathit{\tau }\right)}\mathrm{d}a\mathrm{d}\mathit{\tau }\\ {\displaystyle }{\mathbf{B}}_r=& {\displaystyle }\underset{\mathrm{0}}{\overset{r}{\int }}\underset{\mathit{\tau }}{\overset{r}{\int }}\underset{\mathit{\tau }}{\overset{r}{\int }}{e}^{\mathbf{K}\left(a-\mathit{\tau }\right)}{\mathbf{LL}}^{\mathbf{T}}{e}^{{\mathbf{K}}^{\mathbf{T}}\left(b-\mathit{\tau }\right)}\mathrm{d}a\mathrm{d}b\mathrm{d}\mathit{\tau }\\ {\displaystyle }\mathbf{O}=& {\displaystyle }\mathrm{diag}\left(min\left(\mathrm{1},{\mathbf{HMF}}_{\mathit{\delta }}{\mathit{\mu }}_{n-\mathrm{1}}\right)\right)\end{array}$$

2.5 Extension for strong variations in emanation characteristics or discontinuities

The variance σ of the Brownian motion process in Eq. (7) allows a linear dynamical model as shown in Sect. 2.3 and 2.4 to be tuned in a trade-off-like fashion for one of two properties. Predictions are good in times of relatively constant η and low σ, such that the measurement noise is well filtered at the expense of fidelity in response to steep changes in η, or vice versa. In practical situations however, where the humidity can change rapidly and is known to affect the emanation, a period of re-equilibration is induced in the source, before it returns to somewhat stable behavior. Therefore, experimental time series of the radon source spectra typically show properties that are not well addressed by a single such model of linear dynamics due to the described distinct regimes. This kind of situation is strongly related to the tracking of maneuvering targets, a field in which Bayesian recursive estimation is well established. To obtain good filtering performance in the sense of providing relatively smooth estimates with small uncertainty in case of constant η while retaining the ability to quickly react to steep gradients with associated larger estimation uncertainty, one approach suggested by Nadarajah et al. (2012) and Mazor et al. (1998), among others, is to use the interacting multiple model (IMM) recursive estimator. This is commonly used for object tracking and in this work was adopted to refine the procedure outlined in the previous sections. In the IMM, multiple filters operate on the data at the same time, and their output is mixed based on the likelihood of their measurement predictions. In this way a second, discrete, first-order Markov process describing a discrete random variable $s_{t_n}$ that indexes into the several applicable and differently parameterized linear dynamical models is formally introduced and evolves according to some, potentially parameterized, transition matrix Π. As a result, the filtering and smoothing distributions become the compound distributions $p\left(s_n,{\mathit{x}}_n|{\mathit{y}}_{\mathrm{1}:n}\right)$ and $p\left(s_n,{\mathit{x}}_n|{\mathit{y}}_{\mathrm{1}:N}\right)$, respectively, which now also carry probabilistic information regarding the active model index s_n.

These compound distributions can be decomposed into discrete and continuous components which are approximately represented as mixtures of Gaussians. This kind of system is also called a switching linear dynamical system (SLDS). In the SLDS, exact filtering is not computationally feasible (Barber, 2006; Hartikainen and Särkkä, 2012), since the filtering distribution is a mixture of Gaussians whose components are being multiplied by the number of models at each time step resulting in the exponential growth of components. Most approaches for approximation, like the one employed here, replace the resultant mixture at each step of filtering and smoothing with a smaller one, limiting the number of kept components in the Gaussian mixtures to some fixed upper value I. In case of filtering, these types of algorithms are referred to as Gaussian sum filtering (GSF). Figure 1 gives an outline, both in terms of a component-wise (i.e., marginalized only over the mixture component indices i_n) and a marginal view (i.e., marginalized over both the model indices s_n and the component indices i_n) of the applied GSF method for the final SLDS model for the example of two Gaussian components per mixture. The basis of the prediction and correction steps for each combination of active model and prior components, jointly indexed by s_n, s_n−1, and i_n−1, is detailed in Sect. 2.3 and 2.4. By the discretizations carried out in Sect. 2.3 and 2.4, it is thus implicitly assumed that the active model index can only switch after the observation of each spectrum but not during the integration time of each spectrum. The probabilities $p\left(s_n|{\mathit{y}}_{\mathrm{1}:n}\right)$ are also estimated by the algorithm by evaluation of the likelihood of measuring y_n in each component of the prediction step density and marginalization over x_n and associated i_n, respectively. For further details on the GSF method, the reader is directed to Barber (2006). The results of the GSF algorithm are the (unsummed) factors on the right hand side of the following approximation of the decomposed filtering distribution, which consist of said Gaussian mixtures:

$$\begin{array}{}\text{(14)}& p\left({\mathit{x}}_n|{\mathit{y}}_{\mathrm{1}:n}\right)\approx {\sum }_{s_n}p\left(s_n|{\mathit{y}}_{\mathrm{1}:n}\right){\sum }_{i_n}{w}_{s_n,i_n}p\left({\mathit{x}}_n|s_n,i_n,{\mathit{y}}_{\mathrm{1}:n}\right),\end{array}$$

where s_n is the model index at time step t_n, i_n is the component index of the approximating Gaussian mixtures, and $w_{s_n,i_n}$ are the associated Gaussian mixture weights.

Figure 1 Illustration of the iterative computational methods applied in the Gaussian sum filtering (GSF). The filtering distributions (prior and posterior at time steps t_n−1 and t_n, respectively) are displayed as Gaussian mixtures indexed by i_n for each of the models indexed by s_n. The top panel shows a component-wise view, and the bottom panel shows a marginalized (over s_n and i_n) view of the respective compound distributions. The algorithm can be formally divided into the prediction, correction, and approximation steps. The illustrated iteration starts at some specified prior for t₀. For further details, see text and Appendix A2. The output of the algorithm is the approximated decomposed filtering distribution ${\sum }_{s_n}p\left(s_n|{\mathit{y}}_{\mathrm{1}:n}\right){\sum }_{i_n}{w}_{s_n,i_n}p\left({\mathit{x}}_n|s_n,i_n,{\mathit{y}}_{\mathrm{1}:n}\right)$, the input is a prior and a collection of measurements y_1:N with associated time stamps, measurement times, and values for all parameters (see text). The approximate marginal likelihood of the measurement sequence can be computed alongside the filtering.

In this setting, the smoothing distribution is also a compound distribution in which the components of each mixture are multiplied within each smoothing step in a backwards recursive formulation by the number of linear dynamical models. Therefore, smoothing is also only possible approximately, once again, on the basis of approximating each arising Gaussian mixture with a smaller one. One way to approximately obtain smoothed results in the SLDS setting is thus given by the expectation correction (EC) algorithm introduced in Barber (2006) and shown therein to provide state-of-the-art results, both in terms of computational efficiency and accuracy, using the results of the GSF forward pass and performing a backwards recursion through the time series. The EC algorithm requires, analogous to the GSF forward pass, the propagation and correction methods, in the form of the filtering and smoothing steps, for the parameterized linear dynamical system outlined in Sect. 2.3 and 2.4 acting on each Gaussian mixture component. In this case, however, the integrating behavior of the measurements does not lead to any required adjustments, and the applied equations are thus exactly analogous to the classical Rauch–Tung–Striebel smoother that represent a formal reversal of the dynamics. These have also been used in the original presentation of the EC algorithm (Algorithm 5 in Barber, 2006). Figure 2 provides a graphical illustration similar to the filtering method in Fig. 1 of the EC smoothing backward pass. The EC algorithm provided in Barber (2006) was used without further modifications, and for more information on this algorithm, the reader is directed to this work.

While not a main contribution to this paper, for completeness sake and to facilitate possible re-implementation, both the GSF forward pass and the EC backward pass are outlined in pseudo-code in Appendix A2 in the way they have been implemented here.

Figure 2 Illustration of the backward iterative computational steps (1)–(4) in the expectation correction algorithm applied for correcting the results from the Gaussian sum filtering (GSF) outlined in Fig. 1 to the approximate smoothing solution in the switching linear dynamical system. Input to the algorithm are the results from the GSF and associated time stamps and values for all parameters. Output is the decomposed smoothing distribution ${\sum }_{s_n}p\left(s_n|{\mathit{y}}_{\mathrm{1}:N}\right){\sum }_{i_n}{w}_{s_n,i_n}p\left({\mathit{x}}_n|s_n,i_n,{\mathit{y}}_{\mathrm{1}:N}\right)$ as a Gaussian mixture approximation.

To model the two distinct regimes outlined before, two linear dynamical models are used that share their γ values, but for one (cf. Eq. 8), σ is constrained to a small value. Consequently, one of the models corresponds to regimes of changing η and the other to near constant η. This approach allows us via the model index probabilities $p\left(s_n|{\mathit{y}}_{\mathrm{1}:n}\right)$ and $p\left(s_n|{\mathit{y}}_{\mathrm{1}:N}\right)$ obtained from the GSF and the EC algorithms to also probabilistically identify regimes of constant radon emanation within a time series of recorded spectra, even when the retained activity of ²²²Rn is still in a re-equilibration period, as will be illustrated in the experimental results presented later. In that sense, and by construction of the two models, the model index probabilities can be physically interpreted as the probability of the source to currently have stable emanation characteristics. A non-zero σ acting as a regularization term is, however, still needed in the case of the model corresponding to the constant regimes because otherwise the approximations used within both the GSF and the EC algorithm become numerically unstable. The final model contains four unknown parameters, σ,γ, and the components of Π, which is row-wise normalized and whose components are thus parameterized by two parameters. These four parameters are tuned by minimizing the (approximate) negative marginal log likelihood of the measurement series, $-{\sum }_n\ln p\left({\mathit{y}}_n|{\mathit{y}}_{\mathrm{1}:n-\mathrm{1}},\mathit{\sigma },\mathit{\gamma },\mathbf{\Pi }\right)$, that is accessible in the GSF forward pass (see Algorithm 1 in Appendix A2). For the purposes of this initial investigation, we assume the uncertainty arising from uncertain parameters is negligible in light of the counting statistics, the uncertainty encoded in the prior of the initial state, p(x₀), and the components of H.

2.6 Propagation of the detection efficiency uncertainty

In practice, the components of H are not known without uncertainty, and the uncertainty associated with H is the most significant contribution to the combined uncertainty of η in most practical cases. The outlined formalism above is a way to obtain estimates for $p\left({\mathit{x}}_n|{\mathit{y}}_{\mathrm{1}:n},\mathbf{H}\right)$ and $p\left({\mathit{x}}_n|{\mathit{y}}_{\mathrm{1}:N},\mathbf{H}\right)$ or the extended compound distributions $p\left(s_n,{\mathit{x}}_n|{\mathit{y}}_{\mathrm{1}:n},\mathbf{H}\right)$ and $p\left(s_n,{\mathit{x}}_n|{\mathit{y}}_{\mathrm{1}:N},\mathbf{H}\right)$ detailed in the previous section, respectively. For notational simplicity, the following is formulated with respect to the single-model case but applies analogously to the SLDS model. Formally, the inclusion of this systematic measurement uncertainty into the filtering or smoothing result, respectively, is given by Eq. (15), which is not tractable.

$$\begin{array}{}\text{(15)}& p\left({\mathit{x}}_n|{\mathit{y}}_{\mathrm{1}:n}\right)={\int }_{\mathbf{H}}p\left({\mathit{x}}_n|{\mathit{y}}_{\mathrm{1}:n},\mathbf{H}\right)p\left(\mathbf{H}|{\mathit{y}}_{\mathrm{1}:n}\right)\mathrm{d}\mathbf{H}\end{array}$$

In the present case, we assume H to be constant with respect to time and that its distribution is known from previous measurements, independent of the collection y_1:N. Generally, the uncertainty in H associated with, for example, a prior calibration procedure, can reasonably be described by this approximation, if $\frac{\mathrm{d}\mathbf{H}}{\mathrm{d}t}=\mathrm{0}$. This allows rewriting Eq. (15) to obtain Eq. (16); however, the result remains computationally infeasible, since it would require an infinite amount of filtering and smoothing passes.

$$\begin{array}{}\text{(16)}& p\left({\mathit{x}}_n|{\mathit{y}}_{\mathrm{1}:n}\right)={\int }_{\mathbf{H}}p\left({\mathit{x}}_n|{\mathit{y}}_{\mathrm{1}:n},\mathbf{H}\right)p\left(\mathbf{H}\right)\mathrm{d}\mathbf{H}\end{array}$$

Thus, in the setting at hand, $p\left({\mathit{x}}_n|{\mathit{y}}_{\mathrm{1}:n}\right)$, the result of Eq. (16) is given as an infinite mixture of Gaussian mixtures with weights proportional to p(H). Our strategy of choice to approximately include the systematic uncertainty associated with H is to replace the infinite mixture by a finite version, i.e., to compute the filtering and smoothing densities for specific realizations of H and combine them into mixtures weighted with the likelihood of said realization under the density p(H) known from prior calibration measurements. In practice, this is closely related to the selection of sigma points of H within some prior density and using the associated normalized likelihoods as multiplicative weights with the respective filtering and smoothing mixture weights to obtain extended mixtures that approximately reflect the uncertainty associated with H, a method inspired by the unscented Kalman filter (Julier et al., 2000) for approximate filtering in non-linear systems.

2.7 Implementation details

Implementation of the presented algorithms was carried out in Python using the JAX framework (Bradburry et al., 2018), which provides automatic batching and vectorization, just-in-time compilation, and automatic forward- and reverse-mode differentiation. This allows the optimization of the hyperparameters γ, σ, and Π using the ADAM gradient descent minimization (Kingma and Ba, 2014) routine from JAX together with the automatically computed batch gradient of the negative log likelihood, $-{\sum }_n\ln \left({\mathit{y}}_n|{\mathit{y}}_{\mathrm{1}:n-\mathrm{1}},\mathit{\sigma },\mathit{\gamma },\mathbf{\Pi }\right)$, of the GSF forward pass on the entire dataset (Algorithm 1 in Appendix A2). The number of components retained in the arising Gaussian mixtures in GSF and EC is chosen according to available computational resources. Forward and backward step routines for the linear dynamical system given in Eq. (2) (cf. Appendix A1) were implemented as shown in the auxiliary functions for Algorithms 1 and 2 in Appendix A2, respectively. Gaussian mixture reduction as necessary in the GSF forward and EC backward passes (i.e., the approximation steps in the charts in Figs. 1 and 2) was implemented according to a greedy algorithm (based on Kullback–Leibler divergence (Runnalls, 2007). The mean approximation in the EC backward pass, and implementation of the EC backward and GSF forward passes, were directly adopted as presented in Barber (2006). Gaussian mixture weights and model index probabilities are stored and operated on in ln-space (using the log-sum-exp operation for normalization) for improved numerical stability. Routines for the computation of the Matrices F, U, M, C, and B in Eq. (13) were obtained from symbolic computation using SymPy (Meurer et al., 2017), and the symbolic Jordan canonical form of F and were hardcoded afterward.

3 Application to experimental data

Data for this experiment were generated using an electroplated ²²⁶Ra (104.4±0.4) Bq source (Mertes et al., 2020) mounted on an electrically cooled high-purity germanium (HPGe) detector placed inside a 20 m³ climate chamber. γ-ray spectra were recorded over approximately 85 d at intervals of 10 800 s real time. Within the time series, regions of missing measurements are present, i.e., δ values varied. At specific times the relative humidity inside the climate chamber was varied with the intention of inducing changes in the emanation characteristics of the ²²⁶Ra source. Inside the chamber, relative humidity and temperature were measured in close proximity to the source with a SHT-35 sensor (Sensirion).

For each spectrum, counts above 200 keV were summed, a lifetime scaled background count rate (with associated uncertainty that defines R) was subtracted, and the final algorithms described above (as given in Appendix A2) were applied to the resultant time series of count values (i.e., 1D measurement series) as the input data y_1:N, with five Gaussian components per filter in the GSF forward and EC backward passes. Results are illustrated for the filtering in Fig. 1a and for the smoothing in Fig. 1b. Each r was chosen to reflect the real time of each spectrum as provided by the manufacturer's data-acquisition software (Genie 2000, Mirion Technologies), and each δ was computed from the recorded time stamps of acquisition start points. The dead time of the system was accounted for by correcting the derived count values using the dead-time data as provided by the data-acquisition software. The value of H was determined by measurements of a sealed source of similar type and geometry as presented in Mertes et al. (2020). The uncertainty associated with the ²²⁶Ra activity known from previous measurements detailed in Mertes et al. (2020) was encoded in the density of the prior for the state provided to the algorithm. Apart from the ²²⁶Ra activity, a vague Gaussian prior with a diagonal covariance matrix was chosen for p(x₀). Inherently, the model formulation assumes perfect equilibrium between the SLP and ²²²Rn in the source, which is a small approximation on the timescales at hand. The threshold of 200 keV was chosen because ²²⁶Ra emits γ radiation almost entirely below this energy level, such that the spectrum beyond is almost entirely made up of events associated with the SLP in the source and the background radiation. We chose to neglect the information gained from the spectra regarding the ²²⁶Ra activity because it was not found to substantially improve the prior. The summation of spectra is the most straightforward way to utilize the information contained within each spectrum while keeping the dimensionality of y as small as possible. As a result, the prior density for the ²²⁶Ra activity component of the state is retained over the entire dataset, which is why this component of the state is not shown in Fig 3. The last component of the state vectors, $\frac{\mathrm{d}\mathit{\eta }}{\mathrm{d}t}$, is also not shown in Fig. 3 for visual clarity, since it carries no important information and is merely used as a mathematical tool to specify the properties of the stochastic process of η, the main estimation target of this work. The component $\frac{\mathrm{d}\mathit{\eta }}{\mathrm{d}t}$ is therefore not of any practical meaning.

Both the confidence intervals and the median in Fig. 3 were computed from the marginal cumulative density of the Gaussian mixtures using numerical root finding. Additionally, the confidence intervals include a systematic, Gaussian 1 % uncertainty on the specific value of H which was approximately propagated using the derivation in Sect. 2.6 for five distinct realizations of H (${\mathit{\mu }}_{\mathbf{H}},{\mathit{\mu }}_{\mathbf{H}}\pm {\mathit{\sigma }}_{\mathbf{H}},{\mathit{\mu }}_{\mathbf{H}}\pm \mathrm{2}{\mathit{\sigma }}_{\mathbf{H}}$).

Figure 3 Column (a) shows the filtering and column (b) the smoothing results. Input data (red dots in the top-most row; right scale) and estimated residual ²²²Rn activity (left scale) are depicted in the first row. The second row shows the estimated emanating ²²²Rn given by the SLDS deconvolution approach outlined in Sect. 2.5 (black line; deconvolution result) and erroneously calculated from Eq. (4) (red dots). Deviations between the red dots and the filtering and smoothing results in row two are due to the inability of Eq. (4) to account for dynamic behavior. In the third row, the probability for stable regimes under the model, $p\left(s_n|{\mathit{y}}_{\mathrm{1}:n}\right)$ (column a), and $p\left(s_n|{\mathit{y}}_{\mathrm{1}:N}\right)$ (column b), where s_n equals the index of the model for the stable regime (black; i.e., inferred switch points; see Sect. 2.7) and the relative humidity (red; not used as input data) are shown. The fourth row shows the estimated relative standard uncertainty (calculated from the median and the shown confidence interval) of inferred ²²²Rn emanation.

4 Discussion and conclusion

In the present work, we have summarized and explained the limitations of previously available approaches to estimate ²²²Rn emanation through measurements of the short-lived progeny (SLP) retained within the source. As was derived from first principles, those methods to standardize ²²²Rn emanation are limited to sources with stable characteristics within the operational time and, as such, are generally restricted to use in stable environmental conditions. To alleviate this shortcoming, here we developed an alternative approach that directly infers the conditional probability density for the latent ²²²Rn emanation term from spectral time series of the SLP that remains within the source.

During the application of the resultant algorithms to real-world data, the deviations of the steady-state approximation of the previous method (red dots in the second row of Fig. 3) from the estimated true values (black lines in the second row of Fig. 3) become apparent, underpinning the theoretical derivations. In turn, this means that a thorough analysis of data obtained in this way is restricted to models that account for the dynamic nature of the system, which has not previously been reported.

The specific structure of the filtering and smoothing results in the second row of Fig. 3, showing peaked emanation upon increases in humidity, can easily be explained physically as follows, which justifies the results of the applied method. Considering the time series of count data of the SLP within the source (input data; red dots in the first row of Fig. 3), regions can be seen where changes are occurring much more quickly than would be possible based on the well-known radioactive decay kinetics. As was discussed in Sect. 2.2, the time series of counts is theoretically given by a discretely sampled convolved version of the emanation. Hence, peaked emanation must be occurring, such that the observed time series is possible within the theoretically known decay kinetics. Conversely, the drop in humidity and thus emanation at approximately 70 d does not show this behavior, and the observed ingrowth of the counts directly follows the decay kinetics. Apparently, the behavior depends on the direction of change in the emanation characteristics. This is explained by the fact that upon an increase of the effective diffusion coefficient, the source retains more ²²²Rn atoms than the associated equilibrium value, at which point increased outflow can occur for quick re-equilibration. With the realistic assumption of zero back diffusion from the volume into the source, for a change in the opposite direction, the only way to achieve progeny equilibrium is the typical ingrowth of ²²²Rn, which is the exact behavior shown by the deconvolution result but not by the previous method. Upon fresh preparation of an emanation source, at which point no SLP is present in the source but emanation is still considered to be happening, similar count data to the one past 70 d in Fig. 3 may be observed. This behavior was not previously discussed in Linzmaier and Röttger (2013), where the apparent initial drop in emanation (as computed by Eq. 4) resulting from the initial ingrowth of residual ²²²Rn and SLP was seemingly considered to be its true temporal characteristic, implying that Eq. (4) may be applied in such a case.

In constant regimes, results obtained from previously reported methods converge to the values obtained using the deconvolution approach presented here, as illustrated at approximately days 60 to 70 and past 90 d of the data shown in Fig. 3. While the method we present might not seem beneficial in such constant regimes, the recursive Bayesian approach provides a computationally convenient, mathematically coherent, and flexible way to refine the uncertainty upon observation of streaming data (e.g., obtained by continuous operation of spectrometers) also within constant regimes. As such, here we report for the first time the application of a method whereby time series data of an emanation source can be used to derive correct (near) real time values of the emanation, irrespective of the state of the source. Specifically, the use case of this method and our initial motivation is the implementation of surveillance systems for emanation sources based on spectrometric measurements to improve current state-of-the-art realization, and especially the dissemination of the unit Bq m⁻³ for ²²²Rn in the low-level activity concentration regime. With our contribution, and potential extensions thereof, these types of systems will be enabled in the future. Moreover, experimental investigations of the emanation behavior in response to different environmental conditions are drastically improved by our contribution.

To obtain approximate filtering and smoothing algorithms in the context of radioactivity measurements, we extended the well-known computational methods for inference in linear dynamical systems (i.e., Kalman filter and Rauch–Tung–Striebel smoother) with a computationally convenient approximation for the observed Poisson statistics and the integrating behavior of the measurements in Sect. 2.3 and 2.4. In doing so, we demonstrated that integrating measurements results in a Gaussian process with certain covariance with the latent continuous state which retains the convenient closed form of filtering and smoothing through conjugacy in such linear dynamical models. As was shown, the integrating measurements lead to additional additive uncertainty depending on the variance of the Brownian motion which we consider an intuitive result. These results were used in order to construct the final switching linear dynamical system inference algorithms applied during the experiments.

Within the recorded time series, distinct domains were observed in response to the way the humidity in the chamber was modified, which lends itself to the applied switching dynamical system model, differentiating between stable and non-stable regimes. This approach allows smaller uncertainty to be achieved and smoother functional realizations of η in the somewhat stable regimes, but at the same time gives reasonably high uncertainty for the unstable regimes where the deconvolution result relies on only a few data points. A simpler modeling approach relying only on a single linear dynamical system, such as a more classical version of the Kalman filter, cannot produce smooth results for the constant regimes while retaining the ability to react to steep gradients, since both properties are controlled by the variance of the Brownian motion. While all obvious switching points (induced by changes in the relative humidity) within the time series were captured by the algorithm (third row of Fig. 3), the specific autocorrelation we chose to regularize η leads to smearing of the switching point in the backward (i.e., the smoothing) pass. This is indicated by the fact that the model proposes an unstable state of the source even for times before the humidity has undergone the step changes (i.e., in times before a known change in the source properties has occurred). Note that this is not the case for the filtering results. This is can be attributed to the applied symmetric autocorrelation of η (Eq. 7) and its independency from the residual radon activity, and it may be alleviated by asymmetric autocorrelation or non-linear models but at a substantially higher computational cost. Whether our approach translates well to time series of different characteristics (e.g., drift, smooth changes) is as yet unclear and subject to further studies. Nonetheless, the model parameters provide a way to tune the algorithm for different scenarios.

For an approach like the present one to be applicable in metrology, uncertainty estimates closely related to the guide to the expression of uncertainty in measurement (GUM; Joint Committee for Guides in Metrology, 2008) are needed. At this point, the GUM is restricted to static measurements and first steps are being taken for an extension to dynamic scenarios (e.g., in Eichstädt and Elster, 2012; Elster and Link, 2008; Link and Elster, 2009), where a slightly different formulation for error propagation has been carried out. In the present case, systematic contributions to the uncertainty are dominated by the uncertainty of the measurement mapping through matrix H. We provided a computationally simple approach to approximately propagate this uncertainty across the filtering and smoothing algorithms with arbitrary precision, given that enough computational resources are available and the detection system can reasonably be assumed to be stable in time. Dropping this assumption would require the approximation of the intractable integration in Eq. (15), e.g., through Monte Carlo integration, which was found unnecessary and would have made the algorithm unsuitable for implementation on low-power, portable devices.

Appendix A: Joint density – derivation of discrete forward step for integrating measurement

Assuming that the density $p\left({\mathit{x}}_{n-\mathrm{1}}|{\mathit{y}}_{\mathrm{1}:n-\mathrm{1}}\right)$ is given as $\mathcal{N}\left({\mathit{\mu }}_{n-\mathrm{1}}{\mathbf{\Sigma }}_{n-\mathrm{1}}\right)$, we define the state at the intermediate time point $t_{n-\mathrm{1}}+{\mathit{\delta }}_n$as x_δ for which the following statistics are readily available through the Chapman–Kolmogorov equation (Särkkä and Solin, 2019).

$${\displaystyle }p\left({\mathit{x}}_{\mathit{\delta }}|{\mathit{y}}_{\mathrm{1}:n-\mathrm{1}}\right)\propto \mathcal{N}\left({\mathit{\mu }}_{\mathit{\delta }},{\mathbf{\Sigma }}_{\mathit{\delta }}\right)=\mathcal{N}\left({\mathbf{F}}_{\mathit{\delta }}{\mathit{\mu }}_{n-\mathrm{1}},{\mathbf{F}}_{\mathit{\delta }}{\mathbf{\Sigma }}_{n-\mathrm{1}}{\mathbf{F}}_{\mathit{\delta }}^{\mathbf{T}}+{\mathbf{U}}_{\mathit{\delta }}\right),$$

where F_δ=e^Kδ and ${\mathbf{U}}_{\mathit{\delta }}={\int }_{\mathrm{0}}^{\mathit{\delta }}{\mathrm{e}}^{\mathbf{K}\left(\mathit{\delta }-\mathit{\tau }\right)}{\mathbf{LL}}^{\mathbf{T}}{\mathrm{e}}^{{\mathbf{K}}^{\mathbf{T}}\left(\mathit{\delta }-\mathit{\tau }\right)}\mathrm{d}\mathit{\tau }$, as follows from Eqs. (5) and (11).

By combination of the definition for x(t) in Eq. (10) and the definition of the measurement process for y in Eq. (12), the measurement at t_n is given by the following integration. While the integral in Eq. (12) is an ordinary integral (e.g., in the Riemann sense), the integral in Eq. (10) is a stochastic integral in the Itō sense; i.e., the following double integral is the ordinary integral over an Itō integral, for which we assume that x(t) has continuous sample paths and is square integrable.

$$\begin{array}{ll}{\displaystyle }{\mathit{y}}_{t_n}& {\displaystyle }=\mathbf{H}{\int }_{t_{n-\mathrm{1}}+\mathit{\delta }}^{tn}\mathit{x}\left(\mathit{\tau }\right)\mathrm{d}\mathit{\tau }=\mathbf{H}{\int }_{\mathrm{0}}^r\mathit{x}(t_{n-\mathrm{1}}+\mathit{\delta }+\mathit{\tau })\mathrm{d}\mathit{\tau }\\ {\displaystyle }& {\displaystyle }=\mathbf{H}{\int }_{\mathrm{0}}^r{\mathrm{e}}^{\mathbf{K}s}{\mathit{x}}_{t_{n-\mathrm{1}}+\mathit{\delta }}\mathrm{d}s+\mathbf{H}{\int }_{\mathrm{0}}^r{\int }_{\mathrm{0}}^s{\mathrm{e}}^{\mathbf{K}(s-\mathit{\tau })}\mathbf{L}\mathrm{d}{\mathit{\beta }}_{\mathit{\tau }}\mathrm{d}s\end{array}$$

The triangular domain of the double integral is swapped to find

$${\displaystyle }{\mathit{y}}_{t_n}=\mathbf{H}{\int }_{\mathrm{0}}^r{\mathrm{e}}^{\mathbf{K}s}{\mathit{x}}_{t_{n-\mathrm{1}}+\mathit{\delta }}\mathrm{d}s+\mathbf{H}{\int }_{\mathrm{0}}^r{\int }_{\mathit{\tau }}^r{\mathrm{e}}^{\mathbf{K}(s-\mathit{\tau })}\mathbf{L}\mathrm{d}s\mathrm{d}{\mathit{\beta }}_{\mathit{\tau }},$$

from which the full density $p\left({\mathit{x}}_{n-\mathrm{1}},{\mathit{x}}_{\mathit{\delta }},{\mathit{x}}_n,{\mathit{y}}_n|{\mathit{y}}_{\mathrm{1}:n-\mathrm{1}}\right)$ can be computed using the expectation operator and the definition of the variance, where 〈⋅〉 denotes the expectation operator. Moreover, by definition, $\langle \int \mathrm{d}\mathit{\beta }\rangle =\mathrm{0}$ and thus

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathbf{\Sigma }}_n^{\mathit{y}\mathit{y}}=〈({\mathit{y}}_n-\langle {\mathit{y}}_n\rangle )({\mathit{y}}_n-\langle {\mathit{y}}_n\rangle {)}^{\mathrm{T}}〉\\ {\displaystyle }& {\displaystyle }=〈\mathbf{H}\left({\int }_{\mathrm{0}}^r{\mathrm{e}}^{\mathbf{K}s}\mathrm{d}s{\mathit{x}}_{\mathit{\delta }}-{\int }_{\mathrm{0}}^r{\mathrm{e}}^{\mathbf{K}s}\mathrm{d}s\langle {\mathit{x}}_{\mathit{\delta }}\rangle +{\int }_{\mathrm{0}}^r{\int }_{\mathit{\tau }}^r{\mathrm{e}}^{\mathbf{K}(s-\mathit{\tau })}\mathbf{L}\mathrm{d}s\mathrm{d}{\mathit{\beta }}_{\mathit{\tau }}\right)\\ {\displaystyle }& {\displaystyle }{\left({\int }_{\mathrm{0}}^r{\mathrm{e}}^{\mathbf{K}s}\mathrm{d}s{\mathit{x}}_{\mathit{\delta }}-{\int }_{\mathrm{0}}^r{\mathrm{e}}^{\mathbf{K}s}\mathrm{d}s\langle {\mathit{x}}_{\mathit{\delta }}\rangle +\underset{\mathrm{0}}{\overset{r}{\int }}\underset{\mathit{\tau }}{\overset{r}{\int }}{\mathrm{e}}^{\mathbf{K}(s-\mathit{\tau })}\mathbf{L}\mathrm{d}s\mathrm{d}{\mathit{\beta }}_{\mathit{\tau }}\right)}^{\mathrm{T}}{\mathbf{H}}^{\mathbf{T}}〉.\end{array}$$

Under the assumption of independence of the Brownian motion and x_δ, using Itō isometry and Fubini's theorem, it follows that

$${\mathbf{\Sigma }}_n^{\mathit{y}\mathit{y}}=\langle ({\mathit{y}}_n-\langle {\mathit{y}}_n\rangle )({\mathit{y}}_n-\langle {\mathit{y}}_n\rangle {)}^{\mathrm{T}}\rangle =\mathbf{HM}{\mathbf{\Sigma }}_{\mathit{\delta }}{\mathbf{M}}^{\mathbf{T}}{\mathbf{H}}^{\mathbf{T}}+{\mathbf{HBH}}^{\mathbf{T}},$$

where $\mathbf{M}={\int }_{\mathrm{0}}^r{\mathrm{e}}^{\mathbf{K}s}\mathrm{d}s$ and $\mathbf{B}={\int }_{\mathrm{0}}^r{\int }_{\mathit{\tau }}^r{\int }_{\mathit{\tau }}^r{\mathrm{e}}^{\mathbf{K}\left(a-\mathit{\tau }\right)}{\mathbf{LL}}^{\mathbf{T}}{\mathrm{e}}^{{\mathbf{K}}^{\mathbf{T}}\left(b-\mathit{\tau }\right)}\mathrm{d}a\mathrm{d}b\mathrm{d}\mathit{\tau }$.

Analogously, the cross-covariance between x_n and y_n is given as

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathbf{\Sigma }}_n^{\mathit{x}\mathit{y}}=\langle ({\mathit{x}}_n-\langle {\mathit{x}}_n\rangle )({\mathit{y}}_n-\langle {\mathit{y}}_n\rangle {)}^{\mathrm{T}}\rangle =\\ {\displaystyle }& {\displaystyle }=〈\left({\mathrm{e}}^{\mathbf{K}r}{\mathit{x}}_{\mathit{\delta }}-{\mathrm{e}}^{\mathbf{K}r}\langle {\mathit{x}}_{\mathit{\delta }}\rangle +{\int }_{\mathrm{0}}^r{\mathrm{e}}^{\mathbf{K}\left(r-\mathit{\tau }\right)}\mathbf{L}\mathrm{d}{\mathit{\beta }}_{\mathit{\tau }}\right)\\ {\displaystyle }& {\displaystyle }{\left({\int }_{\mathrm{0}}^r{\mathrm{e}}^{\mathbf{K}s}\mathrm{d}s{\mathit{x}}_{\mathit{\delta }}-{\int }_{\mathrm{0}}^r{\mathrm{e}}^{\mathbf{K}s}ds\langle {\mathit{x}}_{\mathit{\delta }}\rangle +{\int }_{\mathrm{0}}^r{\int }_{\mathit{\tau }}^r{\mathrm{e}}^{\mathbf{K}\left(s-\mathit{\tau }\right)}\mathbf{L}\mathrm{d}s\mathrm{d}{\mathit{\beta }}_{\mathit{\tau }}\right)}^{\mathrm{T}}{\mathbf{H}}^{\mathbf{T}}〉=\\ {\displaystyle }& {\displaystyle }={\mathbf{F}}_{\mathbf{r}}{\mathbf{\Sigma }}_{\mathit{\delta }}{\mathbf{M}}^{\mathbf{T}}{\mathbf{H}}^{\mathbf{T}}+{\mathbf{CH}}^{\mathbf{T}},\end{array}$$

with $\mathbf{C}={\int }_{\mathrm{0}}^r{\int }_{\mathit{\tau }}^r{\mathrm{e}}^{\mathbf{K}\left(r-\mathit{\tau }\right)}{\mathbf{LL}}^{\mathbf{T}}{\mathrm{e}}^{{\mathbf{K}}^{\mathbf{T}}\left(a-\mathit{\tau }\right)}\mathrm{d}a\mathrm{d}\mathit{\tau }$. The result in Eq. (13) is then obtained by marginalizing the full density over x_n−1 and x_δ, which for the Gaussian density means to simply drop the respective columns and rows. Additionally, the forward filtering step is obtained by conditioning the resultant $\mathrm{p}\left({\mathit{x}}_n,{\mathit{y}}_n|{\mathit{y}}_{\mathrm{1}:n-\mathrm{1}}\right)$, Eq. (13), onto the observation of y_n using the well-known conditioning formula for the Gaussian distribution. An implementation for such computation is given by the function FWD_STEP in Appendix A2.

Appendix B: Applied algorithms in pseudo-code