Evaluation of the Long-Term Stability of Metrology Instruments

Abstract

This chapter aims to emphasize the issue of the long-term stability of instruments used in metrology. This issue is a concern mentioned in the IEC/ISO17025:2017 standard and the JCGM100:2008 guide. Control charts are mentioned in these key documents as tools to assess whether a measurement process is under statistical control or not. Control charts (Shewhart charts, CUSUM chart, EWMA chart) are introduced and tested with simulated and real datasets from metrology instruments that operate at the ionizing department of the BIPM. The interest and the limits of such statistical analysis are discussed. They take their basis in a measurement model composed of Gaussian white noise. Although a measurement monitored over a relatively short period may be consistent with this model, it has been observed that the autocorrelation of the measurement data acquired over a long period limits the relevance of control charts. In this case, time series analysis seems more appropriate than conventional control charts. As an illustration, an optimal Bayesian smoother is introduced to demonstrate how to deconvolve the low-frequency random noise and refine the evaluation of uncertainty according to the measurement model for long-term measurement.

1. Introduction

Long-term reproducibility in experimental sciences is a critical metrological topic recently discussed in [1]. In laboratories, the quality system is the frame in which such an issue is controlled. It is notably required in the ISO/IEC 17025 [2] that laboratories must ensure the validity of their results. For this, monitoring of the state of operating of a measurement system must be set up by periodic checks and statistical analysis of the obtained process control value.

Techniques to assess the reproducibility of a process is a matter of concern since the industrial era. An important step has been taken in 1924 with the invention of control charts by Walter Andrew Shewhart (Bell Telephone). Control charts are efficient tools to control if a process is under statistical control and is now wildly used in many domains and recognized by the Project Management Institute as one of the seven basic tools of quality [3]. Moreover, control charts are explicitly proposed in ISO/IEC 17025 [2] and the JCGM100:2008 Guide [4] as tools to implement to analyze periodic checks using control standards.

The main advantages of control charts are their simplicity of implementation and their generic use, independent of the domains. A control chart displays a set of points and control limits. This representation makes it easy to check the number of points exceeding the control limits. Therefore, control charts can be analyzed by anybody without any specific skills in statistics.

Four typical control charts are introduced in this document: the original Shewhart charts (x-chart and R-chart), the CUSUM chart, and the EWMA chart. The manners that they are introduced below should be considered as examples that can slightly vary from one reference to another [5, 6]. These techniques are applied below to simulated data and real data from metrology instruments used in ionizing radiation metrology. Then, the interest of time series analysis is highlighted for the monitoring of measures over a long period using the same datasets.

2. Introduction to control charts

2.1 The x-chart

Periodical control using a check standard consists of periodically reproduce the measurement $x_i$ at the time $t_i$ and forms the measurement set $\left(x_1,x_2,\dots ,x_i,\dots ,x_N\right)$ where $x_N$ is the last recorded point. In the x-chart, a central line is drawn by the moving empirical average ${\overline{x}}_N$ over a period of 31 samples such as

$$\left\{\begin{array}{l}{\overline{x}}_N=\frac{{\sum }_{i=1}^N{x}_i}{N},\text{for}N\le 31\\ {\overline{x}}_N=\frac{{\sum }_{i=N-31}^N{x}_i}{31},\text{for}N>31.\end{array}\right.$$E1

The control limits (LCL and UCL) are given by the multiple of the standard deviations:

$$\begin{array}{c}\left\{\begin{array}{l}s\left(x_N\right)=\sqrt{\frac{{\sum }_{i=1}^N{\left(x_i-{\overline{x}}_N\right)}^2}{N-1}},\text{for}N\le 31,\\ s\left(x_N\right)=\sqrt{\frac{{\sum }_{i=N-31}^N{\left(x_i-{\overline{x}}_N\right)}^2}{30}},\text{for}N>31,\end{array}\right.\end{array}$$E2

such as

$$\begin{array}{c}\left\{\begin{array}{l}\mathrm{LCL}={\overline{\mathrm{x}}}_{\mathrm{N}}-\mathrm{ks}\left({\mathrm{x}}_{\mathrm{N}}\right)\\ \mathrm{UCL}={\overline{\mathrm{x}}}_{\mathrm{N}}+\mathrm{ks}\left({\mathrm{x}}_{\mathrm{N}}\right),\end{array}\right.\end{array}$$E3

where the multiplication factor $k$ is set equal to 2 for the first control level and 3 for the second control level.

Failure criteria are defined according to a white Gaussian model that will be introduced later. It is important to note that the relevance of the analysis is low when the number of measures $N<10$, relevant for $10\le N\le 31,$ and good for $N\ge 31$. Some criteria can be founded in the ISO 7870 standard [5]. As an example, the following rules are proposed where a measurement process is not considered under statistical control when:

• more than 0.35% of the points are above or below the control limits with the multiplication factor $k=3,$

• more than 20% of the points are above or below the control limits with the multiplication factor $k=2$,

• more than 9 points are in a row one side of ${\overline{x}}_N$,

• more than 6 points are in a row increasing or decreasing,

• more than 14 points are in a row alternating between increase and decrease.

2.2 The R-chart

In the R-chart, the plotted points $R_i$ are the derivatives of the measurements:

$$\begin{array}{c}{R}_i=x_i-x_{i-1}\text{for}\mathrm{all}i\in ⟦2,N⟧.\end{array}$$E4

The central line and two levels of control limits are obtained in the same manner as for the x-chart by calculating the moving empirical average and standards deviations of the values $R_i$. Also, the same failure criteria can be applied.

2.3 The CUSUM chart

The cumulative of deviations $C_i$ is obtained by the following recursive formula:

$$\begin{array}{c}{C}_i=C_{i-1}+x_i-{\overline{x}}_N\end{array}$$E5

In this chart, no central line is displayed but two control limits $V_i$ are defined by the so-called V-mask:

$$\begin{array}{c}{V}_i=C_N\pm \mathit{ks}\left(x_i\right)\left(d+N-i\right)\end{array}$$E6

where $k=0.1$ and $d=147.22$.

Any points out of the V-mask mean that a significant drift of the process control value is detected.

2.4 The EWMA chart

The EWMA chart is obtained by plotting the exponentially smoothed data points $z_i$ by using the following recursive formula:

$$\begin{array}{c}{z}_i=\lambda x_i+\left(1-\lambda \right)z_{i-1}\end{array}$$E7

And where $\lambda$ is a smoothing parameter.

The control limits $\mathit{CL}$ are defined by

$$\begin{array}{c}{\mathit{CL}}_i={\overline{x}}_i\pm 3s\left(x_i\right)\sqrt{\frac{\lambda }{2-\lambda }\left(1-{\left(1-\lambda \right)}^{2i}\right)}\end{array}$$E8

With the smoothing parameter $\lambda =2/N$.

Any point out of these control limits means that a significant drift of the process control value is detected.

It must be noted that these two families of control charts, the Shewhart charts on the one hand and the CUSUM and EWMA charts on the other hand, are complementary. Indeed, the Shewhart charts aim to detect only large shifts above $1.5s\left(x_i\right)$ (detection of outliers), whereas the CUSUM and EWMA charts aim to detect only small shifts below $1.5s\left(x_i\right)$ and more generally drifts in the process. In other words, Shewhart charts perform an analysis in the high-frequency domain, whereas the CUSUM and EWMA charts are more suitable to analyze in the low-frequency domain.

3. Test of control charts with simulated data

3.1 The Gaussian white noise model

$N=100$ normally distributed data points $x_i$ are simulated using the random generator from the numpy.py library.

The model of the simulated process is the following:

$$\begin{array}{c}\left\{\begin{array}{c}{x}_i={\mu }_i+{\epsilon }_i\\ {\epsilon }_i\sim \mathcal{N}\left(0,{\sigma }_i^2\right)\end{array}\right.\end{array}$$E9

with ${\mu }_i$: the expectation parameter at the time $t_i$, $\mathcal{N}$: the normal distribution, and ${\sigma }_i^2$: the variance parameter of the measurement noise ${\epsilon }_i$ at the time $t_i$. It is referred below as a Gaussian white noise model.

The reference dataset is defined by setting, for all $i\in ⟦1,N⟧$, the parameters ${\mu }_i=0$ and ${\sigma }_i^2=1$. The control charts are presented below, and the tests conclude—as expected—that the simulated process is under statistical control (Figures 1–4).

Figure 1.

x-chart of the process simulated form the Gaussian white noise model.

Figure 2.

R-chart of the process simulated form the Gaussian white noise model.

Figure 3.

CUSUM chart of the process simulated form the Gaussian white noise model.

Figure 4.

EWMA chart of the process simulated form the Gaussian white noise model.

3.2 The Gaussian white noise model with a deterministic drift

Now, a slight positive drift is considered in the following model referred to as a deterministic drift model in a Gaussian white noise such as the expectation parameter takes values such as:

$$\begin{array}{c}\left\{\begin{array}{ll}{\mu }_i=0,& \forall i\in ⟦1,\frac{N}{2}⟧\\ {\mu }_i={\mu }_{i-1}+0.3,& \forall i\in ⟦\frac{N}{2}+1,N⟧.\end{array}\right.\end{array}$$E10

This leads to a detection of this trend in the x-chart because more than nine points are in a row on one side of the central line. Also, the CUSUM chart and the EWMA chart detect the trend with some points out of their control limits (Figures 5–8).

Figure 5.

x-chart of the process simulated form the Gaussian white noise model with a slight drift.

Figure 6.

R-chart of the process simulated form the Gaussian white noise model with a slight drift.

Figure 7.

CUSUM chart of the process simulated form the Gaussian white noise model with a slight drift.

Figure 8.

EWMA chart of the process simulated form the Gaussian white noise model with a slight drift.

3.3 The autocorrelated model

The third model studied here is the autocorrelated model (or the colored noise model) where a low-frequency Brownian noise is introduced in the following autocorrelated model:

$$\begin{array}{c}\left\{\begin{array}{c}{x}_i=x_{i-1}+{\dot{x}}_i+{\epsilon }_i\\ {\dot{x}}_i={\dot{x}}_{i-1}+{\delta }_i\end{array}\right.\end{array}$$E11

where the measurement noise is ${\epsilon }_i\sim \mathcal{N}\left(0,{\sigma }_i^2\right)$ (idem as in Eq. (9)) and a process noise ${\delta }_i\sim \mathcal{N}\left(0,{\tau }_i^2\right)$ is added with a variance parameter ${\tau }_i^2$. In the example below, ${\tau }_i=0.3$ while ${\sigma }_i$ is still equal to 1.

Again, this type of deviation from the Gaussian white noise model is detected by the value chart with more than nine points in a row on one side of the central line and by CUSUM and EWMA charts (Figures 9–12).

Figure 9.

x-chart of the process simulated form the autocorrelated model.

Figure 10.

Range chart of the process simulated form the autocorrelated model.

Figure 11.

CUSUM chart of the process simulated form the autocorrelated model.

Figure 12.

EWMA chart of the process simulated form the autocorrelated model.

These three examples underline the fact that control charts are—by definition—build to verify whether a dataset fits with a Gaussian white noise model or not. Any deviation from this model is efficiently revealed by control charts. An operator could therefore conclude by the fact that the process is not under statistical control. This conclusion is good in the case of deterministic trends that are that a metrology system seeks to avoid maintaining the references that it produces. However, it has been demonstrated that the autocorrelated noise could also lead to an alert by control charts. On the contrary to the deterministic drift, the latter is a natural random process that appears in every system, which is observed on enough long-term periods. This particular limit of control charts was already pointed out by many authors [7, 8, 9, 10]. This issue is an emphasis in the real cases presented below.

4. Test of control charts with real data

Two data are here considered. The first one comes from the periodic control of a recent metrology system that operates at the BIPM for radionuclide metrology and is called the ESIR (extended SIR) [11, 12]. The second example comes from the periodic control of a metrology service, the SIR, which operates at the BIPM for radionuclide metrology since 1976 [13].

4.1 Test with 1-year period data

The ESIR is under development, and the periodical control using check standards has been set up for 1 year. Figures 13–16 show the analysis of the measured data by control charts. The monitored value is a comparison indicator obtained by using a check source of ¹⁴C. The analysis from the four-control charts leads to the conclusion that the process is under statistical control. In other words, the measurement model complies with a Gaussian white noise model. At least, it is true for this short period of observation. The activity estimation is of the ¹⁴C check source is stable and equal to $24.3373\pm 0.0184$ A.U. There is no unit associated with the measurand because the aim of the ESIR is only to deliver a measure proportional to the activity but is very reproducible during the time. The relative standard uncertainty associated with this measurement is about 0.08%.

Figure 13.

x-chart of the periodic control of the ESIR metrology system.

Figure 14.

R-chart of the periodic control of the ESIR metrology system.

Figure 15.

CUSUM chart of the periodic control of the ESIR metrology system.

Figure 16.

EWMA chart of the periodic control of the ESIR metrology system.

4.2 Test with 20-year period data

The SIR produces international references for decades. The ²²⁶Ra source is used as a control standard to monitor the reproducibility of the measured current. In this case, the control charts are applied to a 20-year period of data acquisition. They have detected that the process is not under statistical control. However, these low-frequency fluctuations could certainly be a random colored noise and not a deterministic drift. This phenomenon appears when measurement values are observed over a very long period such as here. The only conclusion that can be drawn from control charts is that the data do not fit with the Gaussian white noise model and that a type A evaluation of the uncertainty associated with the measurement values (from t₁ to t_N) cannot be applied. Nevertheless, it is no need to refine the measurement model for the SIR measurement because the quantity delivered by the system is a ratio of a punctual current measurement (at a time t_i) from which this low-frequency component revealed here is suppressed by the method (Figure 17–20).

Figure 17.

x-chart of the periodic control of the SIR metrology system.

Figure 18.

R-chart of the periodic control of the SIR metrology system.

Figure 19.

CUSUM chart of the periodic control of the SIR metrology system.

Figure 20.

EWMA chart of the periodic control of the SIR metrology system.

We have underlined, through these simulated and real datasets, the limitation of control charts to deal with the evaluation of the long-term stability of metrology instruments. In this case, time series analysis could be used as more appropriate tool to address this objective.

An example of time series analysis is introduced below.

5. Introduction to time series analysis applied to a long-term measurement model

It has been seen that the autocorrelated model introduced in Eq. (10) is appropriate to represent measurement data when observed over a long-term period. Eq. (10) is here rewritten in a state-space model where the measurement is decomposed by the following:

• the fluctuation of the expectation value $\left(x_1,\dots ,x_i,\dots ,x_N\right)$ with the low-frequency random noise at measurement times $\left(t_1,\dots ,t_i,\dots ,t_N\right)$,

• and the measurements gathered in the vector $\left(y_1,\dots ,y_i,\dots ,y_N\right)$ where the Gaussian white noise is finally added to $\left(x_1,\dots ,x_i,\dots ,x_N\right)$.

The model becomes,

$$\begin{array}{c}\left\{\begin{array}{l}{x}_i=x_{i-1}+{\dot{x}}_{i-1}\\ {\dot{x}}_i={\dot{x}}_{i-1}+{\delta }_i\\ y_i=x_i+{\epsilon }_i\end{array}\right.\end{array}$$E12

The drift velocity ${\dot{x}}_i=x_i-x_{i-1}$ is introduced such that the dynamic Gaussian noise ${\delta }_i\sim \mathcal{N}\left(0,{\tau }^2\right)$ describes a random walk parametrized by a variance ${\tau }^2$. And the measurement noise ${\epsilon }_i\sim \mathcal{N}\left(0,{\sigma }^2\right)$ is Gaussian white noise with a variance ${\sigma }^2$.

The idea behind time series analysis is to deconvolve the low-frequency fluctuation of the expectation values $\left(x_1,\dots ,x_i,\dots ,x_N\right)$ from the observations $\left(y_1,\dots ,y_i,\dots ,y_N\right)$.

At first, the measurement variance ${\sigma }^2$ is evaluated by using the median of all the punctual type A evaluation of the standard deviation $s\left(y_i\right)$ which are (by definition) isolated from the dynamic component.

$$\begin{array}{c}{\sigma }^2=\text{median}\left(s^2\left(y_i\right)\right)\end{array}$$E13

with,

$$\begin{array}{c}{s}^2\left(y_i\right)=\frac{{\sum }_{j=1}^M{\left(y_i^{\left(j\right)}-M^{-1}{\sum }_{l=1}^M{y}_i^{\left(l\right)}\right)}^2}{M-1}.\end{array}$$E14

$\left(y_{N+1}^{\left(1\right)},\dots ,y_{N+1}^{\left(j\right)},\dots ,y_{N+1}^{\left(M\right)}\right)$ is the set of measurements obtained at the time $t_i$ during a period that allowing type A evaluation according to the Gaussian white noise model. Here, control charts can be used to verify this hypothesis.

Eq. (12) is transformed in the following matrix form:

$$\begin{array}{c}\left\{\begin{array}{l}{X}_i=F{X}_{i-1}+W_i\\ y_i=H{X}_i+{\epsilon }_i\end{array}\right.\end{array}$$E15

with the “state vector” $X_i=\left(\genfrac{}{}{0pt}{}{x_i}{{\dot{x}}_i}\right)$, the “transition matrix” $F=\left(\begin{array}{cc}1& 1\\ 0& 1\end{array}\right)$, the” dynamic noise matrix” $W_k=\left(\genfrac{}{}{0pt}{}{0}{{\delta }_i}\right),$ and the “observation matrix” $H=\left(\begin{array}{cc}1& 0\end{array}\right)$.

5.1 The Kalman filter

Given the Gaussianity of the noises (both dynamic noise ${\delta }_i$ and measurement noise ${\epsilon }_i$) of this linear model (in the sense of time series), the Kalman filter (KF) provides an optimal estimation of the states $X_i$ given the current and past measurements $y_{1:i}$. The optimality of KF to solve Eq. (15) is established through both least-square interpretation and the probabilistic Bayesian analysis [14, 15]. Hence, KF is a recursive Bayesian filter where the states $X_i$ and the measurements $y_i$ are Markov processes implying

$$\begin{array}{c}\left\{\begin{array}{c}p\left(X_i|X_{1:i-1}\right)=p\left(X_i|X_{i-1}\right)\\ p\left(y_i|X_{1:i}\right)=p\left(y_i|X_i\right).\end{array}\right.\end{array}$$E16

KF evaluates, at each time step $i$, the parameters of the Gaussian posterior distribution $p\left(X_i|y_{1:i}\right)$ in two steps:

• The first step consists of predicting the state $X_i$ given past measurement data $y_{1:i-1}$. The posterior distribution $p\left(X_{i-1}|y_{1:i-1}\right)$ at $t_{i-1}$ is used to build a prior distribution at the time $t_i$ according to the linear model (see Eq. (15)) from $t_{i-1}$ to $t_i$ such as

$$\begin{array}{c}p\left(X_i|y_{i-1}\right)=\int p\left(X_i|X_{i-1}\right)p\left(X_{i-1}|y_{i-1}\right)\mathrm{d}{X}_{i-1}.\end{array}$$E17

• The second step aims to update the posterior distribution at $t_i$ given the predicted prior $p\left(X_i|y_{1:i-1}\right)$, and the likelihood $p\left(y_i|X_i\right)$ of the measurement $y_i$ given the estimate $X_i$.

$$\begin{array}{c}p\left(X_i|y_i\right)\propto p\left(X_i|y_{i-1}\right)p\left(y_i|X_i\right).\end{array}$$E18

The posterior state estimate ${\widehat{X}}_{i\mid i}$ and the posterior of the covariance $P_{i\mid i}$ of the posterior distribution $p\left(X_i|y_{1:i}\right)$ are evaluated through the Kalman filtering process described below.

From $i=1$ to $i=N$, KF predicts the state $X_i$ and the measurement $y_i$ at the time step $i$ given the previous state at $i-1$:

$$\begin{array}{c}\left\{\begin{array}{c}{\widehat{X}}_{i\mid i-1}=F{\widehat{X}}_{i-1\mid i-1}\\ {\widehat{y}}_{i\mid i-1}=H{\widehat{X}}_{i\mid i-1}.\end{array}\right.\end{array}$$E19

With the state and measurement covariances, $P_i$ and $S_i$ such that:

$$\begin{array}{c}\left\{\begin{array}{c}{P}_{i\mid i-1}=F{P}_{i-1\mid i-1}{F}^{\mathrm{T}}+{\tau }^2\\ S_i=H{P}_{i\mid i-1}{H}^{\mathrm{T}}+{\sigma }^2.\end{array}\right.\end{array}$$E20

Then, the pre-fit residual (innovation) $e_i$ of the estimation and the Kalman gain $K_i$ are calculated by

$$\begin{array}{c}\left\{\begin{array}{c}{e}_i=y_i-{\widehat{y}}_{i\mid i-1}\\ K_i=P_{i\mid i-1}{H}^{\mathrm{T}}{S}_i^{-1}.\end{array}\right.\end{array}$$E21

Finally, the updated state estimation is

$$\begin{array}{c}\left\{\begin{array}{c}{\widehat{X}}_{i\mid i}={\widehat{X}}_{i\mid i-1}+K_i{e}_i\\ P_{i\mid i}=\left(I-K_{i}H\right)P_{i\mid i-1},\end{array}\right.\end{array}$$E22

with $I$: the identity matrix.

5.2 The Rauch-Tung-Striebel smoother

The KF estimates the states ${\widehat{X}}_{i\mid 1:i}$ given present and past observations from $t_1$ to $t_i$. To refine the estimation, a KF backward pass (smoothing pass) processes the filtered estimates data from $t_N$ to $t_1$. This procedure is called the Rauch-Tung-Striebel smoother (RST) [16, 17]. Symmetrically to the forward KF optimal filter, the optimal smoothing also consists in:

• a prediction step where the state $X_{i+1}$ is estimated given measurement data $y_{1:i}$ and the forecasted state backward transition according to the linear model from $t_{i+1}$ to $t_i$.

$$\begin{array}{c}p\left(X_{i+1}|y_{1:i}\right)=\int p\left(X_{i+1}|X_i\right)p\left(X_i|y_{1:i}\right)\mathrm{d}{X}_i\end{array}$$E23

• The second step aims to update the posterior distribution at $t_i$ given all the data $p\left(X_i|y_{1:N}\right)$.

$$\begin{array}{c}p\left(X_i|y_{1:N}\right)=p\left(X_i|y_{1:i}\right)\int \frac{p\left(X_{i+1}|X_i\right)p\left(X_{i+1}|y_{1:N}\right)}{p\left(X_{i+1}|y_{1:i}\right)}\mathrm{d}{X}_{i+1}.\end{array}$$E24

The RTS smoother calculates the Maximum A Posteriori (MAP) estimate ${\widehat{X}}_{i\mid 1:N}$ of the posterior probability $p\left(X_i|y_{1:N}\right)$ given all the observations (past, present, and future) from $t_1$ to $t_N$.

From $i=N$ to $i=1$, the residual $e_i^{\prime }$ between the filtered state ${\widehat{x}}_{i+1\mid i}$ and the prediction ${\widehat{x}}_{i+1\mid 1:N}$ is used to calculate a new Kalman gain $G_i$ such as

$$\begin{array}{c}\left\{\begin{array}{l}{e}_i^{\prime }={\widehat{X}}_{i+1\mid 1:N}-{\widehat{X}}_{i+1\mid i}\\ G_i=P_{i\mid i+1}{F}^{\mathrm{T}}{P}_{i+1\mid i}^{-1}.\end{array}\right.\end{array}$$E25

The optimal state estimation is finally given by

$$\begin{array}{c}\left\{\begin{array}{l}{\widehat{X}}_{i\mid 1:N}={\widehat{X}}_{i\mid i}-G_i{e}_i^{\prime }\\ P_{i\mid 1:N}=P_{i\mid i}+G_i\left(P_{i+1\mid 1:N}-P_{i+1\mid i}\right)G_i^{\mathrm{T}}.\end{array}\right.\end{array}$$E26

5.3 The auto-tuning of the RTS smoother

As the measurement noise variance ${\sigma }^2$ is evaluated, by definition, through measurement repeated during a period free from long-term fluctuation (Eq. (17)), the variance of the dynamic noise ${\tau }^2$ is the remaining unknown parameter of the model. The auto-tuning procedure proposed here to evaluate $\tau$ starts with a loss function $L$ defined by

$$\begin{array}{c}L\left(\tau \right)=\sum _{i=2}^N{e}_i^2+{\left({\widehat{x}}_{i\mid 1:N}-{\widehat{x}}_{i-1\mid 1:N}\right)}^2.\end{array}$$E27

The value $\widehat{\tau }$ minimizing $L\left(\tau \right)$ aims to minimize the residual between the measurements $y_i$ and the estimations of the measurements ${\widehat{y}}_{i\mid i-1}$ as well as the distance between two consecutive filtered estimations ${\widehat{x}}_{i\mid 1:N}$ and ${\widehat{x}}_{i-1\mid 1:N}$. On the one hand, if this last smoothing constraint is not considered, $\tau$ will diverge toward high values describing a totally unstable system. On the other hand, if the prediction error of the model is not considered, $\tau$ will be canceled by describing a stable system with constant ${\widehat{x}}_{i\mid i}$ values. So, the loss function of Eq. (27) permits to produce a model of measurement of a controlled process containing a low-frequency dynamic noise. The variance of the low-frequency dynamic noise is determined through an iterative procedure minimizing $L\left(\tau \right)$ such that,

$$\begin{array}{c}\widehat{\tau }=\underset{\tau }{\text{argmin}}\left(L\left(\tau \right)\right).\end{array}$$E28

The auto-tune RTS smoother provides filtered values ${\widehat{x}}_{i\mid 1:N}$ where the measurement noise, ${\epsilon }_i,$is removed given the knowledge of the variances of the measurement noise ${\sigma }^2$ (see Eqs. (13) and (14)) and the dynamic noise ${\widehat{\tau }}^2$. So, the global hidden uncertainty component from the long-term process random noise can be evaluated from the deconvolved data ${\widehat{x}}_{i\mid 1:N}$ such that,

$$s\left({\widehat{x}}_{i\mid 1:N}\right)={\left(N-1\right)}^{-1}\sum _{i=1}^N{\left({\widehat{x}}_{i\mid 1:N}-N^{-1}\sum _{j=1}^N{\widehat{x}}_{i\mid 1:N}\right)}^2.$$E29

Thanks to the intrinsic unfolding of the measurement noise (Gaussian white) and the process noise (Gaussian red) defined in the dynamic model (Eq. 12), the standard uncertainty $u\left(x_{N+1}\right)$ of a new measurement at $t_{N+1}$ can be obtained by the quadratic summation of the type A uncertainty evaluation $s\left(x_{i\left(N+1\right)}\right)$ obtained under repeated conditions $\left(y_{N+1}^{\left(1\right)},\dots ,y_{N+1}^{\left(j\right)},\dots ,y_{N+1}^{\left(M\right)}\right)$ and the type B uncertainty evaluation evaluated using all historical data $\left(y_1,\dots ,y_i,\dots ,y_N\right)$ through the above-described auto-tunned RTS smoother.

Therefore, the standard uncertainty $u^2\left(x_{N+1}\right)$ of the next measurement at $t_{N+1}$ should be inflated considering the impact of the dynamic noise such that,

$$u^2\left(x_{T+1}\right)=s^2\left(x_{i\left(T+1\right)}\right)+P_{N+1\mid 1:N}\left(1\right)$$E30

With,

$$\begin{array}{c}{P}_{N+1\mid 1:N}=F{P}_{N\mid 1:N}{F}^{\mathrm{T}}+{\tau }^2\end{array}$$E31

This auto-tuning of a RTS smoother applied to measurement data recorded over a long-term period fully address the recommendation of the JCGM100:2008 guide [4]:

“Because the mathematical model may be incomplete, all relevant quantities should be varied to the fullest practicable extent so that the evaluation of uncertainty can be based as much as possible on observed data. Whenever feasible, the use of empirical models of the measurement founded on long-term quantitative data, and the use of check standards and control charts that can indicate if a measurement is under statistical control, should be part of the effort to obtain reliable evaluations of uncertainty. The mathematical model should always be revised when the observed data, including the result of independent determinations of the same measurand, demonstrate that the model is incomplete. A well-designed experiment can greatly facilitate reliable evaluations of uncertainty and is an important part of the art of measurement.”

6. Application of time series analysis on simulated and real measurement data

The time series analysis through the auto-tuned RST smother is presented in Figure 21 for the simulated dataset from Section 3.3 obtained with the autocorrelated model. Black points are the measurement points $y_i$ with error bars representing the standard deviation $s\left(y_i\right)$. Green points are the estimations of the low-frequency fluctuations (red noise) of the measurement ${\widehat{x}}_i$ with error bars representing the standard deviation $s\left({\widehat{x}}_i\right)$. The red line is the global expectation value where the thickness of the line represents +/− 1 standard deviation.

Figure 21.

Time series analysis of the measurement data simulated form the autocorrelated model.

It can be seen the good deconvolution of the low-frequency random noise shown in green points. This leads to a standard uncertainty from this autocorrelated noise, $s\left({\widehat{x}}_{i\mid 1:N}\right)=1.31.$ This value is significant when compared to type A evaluation of individual measurement points ${\sigma }^2=1.74$.

Also, the time series analysis through the auto-tuned RST smother is presented in Figure 22 when applied to the periodic control of the SIR metrology system. Again, black points are the measurement points $y_i$ with error bars representing the standard deviation $s\left(y_i\right)$. Green points are the estimations of the low-frequency fluctuations (red noise) of the measurement ${\widehat{x}}_i$ with error bars representing the standard deviation $s\left({\widehat{x}}_i\right)$. The red line is the global expectation value where the thickness of the line represents +/− 1 standard deviation.

Figure 22.

Time series analysis of the periodic control of the SIR metrology system.

The analysis using control charts (see Section 4.2) has led to conclude that the system is not under statistical control although autocorrelated noise is expected for such a long-term period of monitoring. The optimal Bayesian smoother permits to properly deconvolve this component (see green points) and to evaluate the standard uncertainty associated with it: $s\left({\widehat{x}}_{i\mid 1:N}\right)=0.0183$. This value is significant when compared to type A evaluation of individual measurement points ${\sigma }^2=0.025$.

7. Conclusion

It has been emphasized through simulated and experimental data that control charts (Shewhart, CUSUM, EWMA charts) are convenient tools to assess whether a measurement system is under statistical control or not. They test the compliance of the measured dataset with a Gaussian white noise model. It can notably be used to determine the possibility to apply a type A evaluation of uncertainty on the measurement data.

However, care should be taken when monitoring over long periods of time. A low-frequency random noise will inevitably appear. Despite its stochastic nature, the latter will induce detection by control charts. In this case, it is not appropriate to conclude by a loss of control of the system but rather to re-evaluate the measurement model to include this long-term component. To that end, time series analysis could be implemented to better deal with this particular case of measurement data.

An auto-tuned Rauch-Tung-Striebel smoother, based on the optimal Bayesian filter called Kalman filter, is introduced to illustrate how time series analysis can help to deconvolve the components of the long-term measurement model. The low-frequency random fluctuation (aka red noise) can be estimated and used to evaluate the standard uncertainty component it induces on the measurement.