A data-driven persistence test for robust (probabilistic) quality control of measured environmental time series: constant value episodes

Robust quality control is a prerequisite and an essential component in any data application. That is especially important for time series of environmental observations such as air quality due to their dynamic and irreversible nature. One of the common issues in these data is constant value episodes (CVEs), where a set of consecutive data values remains constant over a given period. Although CVEs are often considered to be an indicator of sensor failure or other measurement errors and are removed during quality control procedures, there are situations when CVEs reflect natural environmental phenomena, and they should not be removed from the data or analysis. Assessing whether the CVEs are erroneous data or valid observations is a challenge. As there are no formal procedures established for this, their classification is based on subjective judgment and is therefore uncertain and irreproducible. This paper presents a novel test procedure, i.e., constant value test, to estimate the probability of CVEs being valid data. The theoretical foundation of this test is based on statistical characteristics and probability theory and takes into account the numerical precision of the data values. The test is a data-driven (parametric) approach, which makes it usable for time series analysis in different environmental research domains, as long as serial dependency is given and the data distribution is not too different from Gaussian. The robustness of the test was demonstrated with sensitivity studies using synthetic data with different distributions. Example applications to measured air temperature and ozone mixing ratio data confirm the versatility of the test.

1 Introduction

Millions of sensors monitor the environment every day, and their data are used in many applications such as trend analysis (Fang et al., 2013; Mills et al., 2016, 2018; Chang et al., 2017; Fleming et al., 2018; Lefohn et al., 2018) and forecasts (Gardner, 1999; Zhang et al., 2012; Debry and Mallet, 2014; Zhou et al., 2019) to provide important information on global challenges such as climate change, air quality, soil degradation, etc. The measurement process can be interpreted as sampling from a true distribution of atmospheric state variables, for example, temperature or air pollutant concentration, at a given location. Each measured value is an estimation of “truth” that has been obtained through a set of data samples (Grant and Leavenworth, 1996). A common feature of many environmental time series is the fact that the true distribution changes with time. This makes such measurements irreproducible.

Measured data can be contaminated by various errors such as systematic, random, non-representative and gross errors (Gandin, 1988; Steinacker et al., 2011). These errors can arise from poor sensor calibration, long-term sensor drift, noise, non-resolvable processes by an observational network, and mistakes during data processing, decoding or transmission. Some of these errors arise from unpredictable natural phenomena such as floods, fire, frost and animal activities (Campbell et al., 2013) that cannot be documented in every detail. Although many efforts are devoted to developing advanced analytical tools and methods, these errors can have deleterious effects on the statistical analyses. For instance, outliers, i.e., values far outside of the norm for a variable or population, can increase the error variance or reduce the power of statistical tests (Osborne and Overbay, 2004). Specifically, constant value episodes (CVEs) can decrease the normality when the assumption of a normal distribution must be satisfied, for example, in linear regression. Thus, even the most sophisticated statistical model can be vulnerable against unknown and potentially erroneous data. If such errors in the data are not identified by applying quality control (QC) procedures, the information obtained from the data will be misleading, and the results from scientific data analyses can be unreliable and biased. Therefore, robust QC procedures are an essential component in the data production chain and a requirement for having a more reliable quantification of trend or other statistical analysis.

Many research initiatives and environmental monitoring programs have thus established standards and guidelines for QC procedures. Most of them rely on visual screening of data, and therefore personal inspection, and on manual elimination of erroneous values based on empirical knowledge and investigator experiences. Several advanced tools such as GCE (Scully-Allison et al., 2018), CoTeDe (Castelão, 2016), and AutoQC (Good et al., 2022) and comprehensive user manuals such as QARTOD (Bushnell et al., 2019) and WMO-AWS (Zahumenský, 2004) have been developed with precise rules to overcome this subjectivity. However, their application is often limited to a few variables or specific data sets, for example, from limited geographic regions with relatively homogenous conditions. This, in turn, can be problematic if one wants to assemble global data sets of various environmental variables. For example, in the Tropospheric Ozone Assessment Report (TOAR), a global database with ground-level ozone measurements at more than 10 000 locations around the world, was built with data from more than 30 different contributors (Schultz et al., 2017). Different QC procedures at these agencies and sites led to increased uncertainty in the assessment. At this scale of data, manual inspection methods are not only error prone but also impractical. It is therefore desirable to develop a more generic, robust and data-driven approach for the QC of environmental monitoring time series.

The focus of this study is to develop a QC test for CVEs as the first element for such data-driven QC. CVEs are a common feature in air quality time series and other environmental data sets. As an example, in a specific 35-year-long ozone time series with hourly sampling, CVEs with a length of 2 occurred 20 313 times. Therefore, about 6.7 % of the data values are CVEs, meaning that such incidents are expected to occur naturally about 16 times per 10 d in the hourly data. The CVEs with a longer length, e.g., 3, 4 and 5, occur 6190, 2887 and 1681 times, respectively, and so the proportion of these incidents are 4.85, 2.26 and 1.31 for 10 d hourly data time series. While they can be detected through a persistence test, a qualified judgment whether such data are erroneous or not is a difficult undertaking. If CVEs are excluded from the data (Horsburgh et al., 2015; Gudmundsson et al., 2018), the results of the analysis, such as model–data comparisons (Bey et al., 2001; Horowitz et al., 2003; Dawson et al., 2008; Emmons et al., 2010; Lamarque et al., 2012; Rasmussen et al., 2012; Tilmes et al., 2012; Im et al., 2015; Schnell et al., 2015; Lyapina et al., 2016; Sofen et al., 2016), can become biased. That can be an issue in (re)analysis products (Inness et al., 2019; Hersbach et al., 2020), where assimilation processes reduce misfits between observations and their modeled values. If the models correctly capture CVE events, excluding the CVEs will lead to a type I error. On the other hand, if CVEs originating from instrument malfunctions are included in the analysis, that will raise type I and type II errors and likely raise unreliable results.

This study presents a new (QC) test procedure, i.e., constant value test (CVT), which estimates the probability of a CVE representing valid data. Data users can select a threshold of an acceptable probability depending on their scientific study or data analysis task. The CVT is entirely data-driven and makes very few assumptions about the properties of the underlying values' distribution and probability density function (Gaussian). Currently, the method is valid for data with a Gaussian frequency distribution. Possible extensions of the method are discussed in the conclusions section. In principle, it is possible to use the technique of statistical simulations to examine how the CVE probabilities change for non-Gaussian distributions. However, this is beyond the scope of this paper. Due to its generality, the test is applicable for a wide variety of environmental variables with a serial dependency (autocorrelation). The article structure is as follows: the method (CVT) is described in Sect. 2. In Sect. 3, the approach is evaluated using synthetic data for demonstration purposes. The results of three real test cases are discussed in Sect. 4. And, finally, conclusions are given in Sect. 5.

2 Methodology

Before describing the proposed method, we briefly summarize some issues with existing methods. In existing QC frameworks, the persistence test is typically defined based on the minimum expected variability, but this requires prior knowledge about the true statistical distribution of the measurements. For example, Zahumenský (2004) has defined that air temperature measurements shall be flagged as “doubtful or suspected value” if the measured variable varies by less than 0.1 K over 60 min. Such a priori assumptions may lead to false data labeling when environmental conditions are exceptionally stable and the true data variability is reduced for some period of time. For instance, temperature variation of 0.1 K can occur in the morning when radiative forcing is small, e.g., on a foggy day in autumn. In measurements of air pollutant concentrations, longer periods of zero values can be found if the measured concentrations are below the instrument detection limit or if chemical conversion leads to a complete removal of a species. For example, ground-level ozone concentrations at urban sites remain zero for several hours if there is a high level of nitrogen oxide emission.

The assessment of CVEs will also have to depend on the numerical precision or resolution (res), which is the number of significant digits with which an observation is recorded (Chapman, 2005). For example, historical measurements of ground-level ozone (Azusa station) in the EPA Air Quality System (AQS) in the 1980s were reported with a resolution of 8 ppb (parts per billion). Another pollutant in the EPA AQS database for which reporting precision has changed over time since 1980 is carbon monoxide at the Fresno station (California state). So, it is not uncommon to find episodes of several hours when all measurements are reported as the same value, and it would be implausible to remove all of them as “erroneous measurements”.

The CVT takes these considerations into account and provides a data-driven approach with very few a priori assumptions. It consists of two main procedures: first, CVEs need to be found and the length of the episodes must be recorded, then, in the second step, the probability of each CVE being a period of valid data with low variability is estimated. While the first procedure can be simply implemented by taking the differences of consecutive values, a possible complication arises if the time series contains missing data or if the data were irregularly sampled. While the software accompanying this paper has a provision to deal with missing data, we ignore the second issue for the purpose of this paper and require that the time series has been sampled at regular intervals. The following method description focuses on the estimation of the likelihood that two or more constant values occur in reality and are thus not necessarily resulting from measurement or data processing errors.

2.1 Statistical background

To describe the joint process of a given time series, we assume such a stochastic process can be represented as a multivariate Gaussian distribution (Tong, 1990; Rencher, 2002). Let $X=(x_{\mathrm{1}},\mathrm{\dots },x_n)$ be a series of random variables; the joint probability density function of a multivariate Gaussian distribution, 𝒩(μΣ), can be written as

$$\begin{array}{}\text{(1)}& f_X\left(x_{\mathrm{1}},\mathrm{\dots },x_n\right)={\displaystyle \frac{\exp \left(-\frac{\mathrm{1}}{\mathrm{2}}(x-\mathit{\mu }{)}^T{\mathrm{\Sigma }}^{-\mathrm{1}}(x-\mathit{\mu })\right)}{\sqrt{(\mathrm{2}\mathit{\pi }{)}^k\mathrm{|}\mathbf{\Sigma }\mathrm{|}}}}.\end{array}$$

Here, μ is an n×1 mean vector and Σ is an n×n positive definite covariance matrix. In the stationary case, without loss of generality, μ can be assumed to be a constant, and Σ can be represented as multiplication of a finite constant variance σ² and a (auto)correlation matrix $\left\{i=\mathrm{1},\mathrm{\dots },n;j=\mathrm{1},\mathrm{\dots },n\right\}$, with ∅(ij)=1 if i=j (diagonal) and $\mathrm{0}\le \mathrm{\varnothing }\left(ij\right)\le \mathrm{1}$ if i≠j (off-diagonal) for a given time series.

Long range approximation of an environmental time series is generally unnecessary and computationally expensive (e.g., Wincek and Reinsel, 1986; Guttorp et al., 1994; Niu, 1996; Fioletov and Shepherd, 2003; Kumar and De Ridder, 2010). Here we use an assumption that an environmental time series is auto-correlated and can be approximated by an autoregressive (AR(1)) process (Tiao et al., 1990; Weatherhead at al., 1998, 2000; Reinsel et al., 2002). The definition of an AR(1) process, the x_i, i.e., data value at time i, can be written as

$$\begin{array}{}\text{(2)}& x_i=\text{const}+\mathrm{\varnothing }{x}_{i-\mathrm{1}}+{\mathit{\epsilon }}_i.\end{array}$$

Here, ε_i is a white noise, and const is an offset. With the assumption of the AR(1) process, the correlation matrix can be approximated by one parameter, ∅, since $\text{Corr}\left(X_i,X_{i-h}\right)={\mathrm{\varnothing }}^{\left|h\right|}$ (the correlation between any two points is only dependent on the time interval, h); thus, the stochastic process can be governed by three parameters, i.e., μ, σ² and ∅.

The general likelihood of an AR(1) process can be approximated using the first-order Markov property as

$$\begin{array}{}\text{(3)}& p\left(x_{\mathrm{1}},\mathrm{\dots },x_n\right)=p\left(x_{\mathrm{1}}\right){\prod }_{k=\mathrm{2}}^{n}p\left(x_k\mathrm{|}{x}_{k-\mathrm{1}}\right),\end{array}$$

where p(x₁) is the density of initial state, which is not critical in this study, because the focus is placed on the probability of a consecutive state that is identical to the previous value, i.e., the second term, and p(x_k | x_k−1) represents the probability distribution of x_k depending only on x_k−1. The above equation is a general form without a distributional assumption. To derive the explicit form for the Gaussian case, we start from a univariate and a bivariate probability density function:

$$\begin{array}{}\text{(4)}& {\displaystyle }f\left(x_{k-\mathrm{1}}\right)={\displaystyle \frac{\mathrm{1}}{\mathit{\sigma }\sqrt{\mathrm{2}\mathit{\pi }}}}\exp \left(-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[{\displaystyle \frac{{\left(x_{k-\mathrm{1}}-\mathit{\mu }\right)}^{\mathrm{2}}}{{\mathit{\sigma }}^{\mathrm{2}}}}\right]\right),\text{(5)}& {\displaystyle }\begin{array}{rl}f\left(x_{k-\mathrm{1}},x_k\right)& ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }{\mathit{\sigma }}^{\mathrm{2}}\sqrt{\mathrm{1}-{\mathrm{\varnothing }}^{\mathrm{2}}}}}\exp (-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}\left(\mathrm{1}-{\mathrm{\varnothing }}^{\mathrm{2}}\right)}}\\ & [{\displaystyle \frac{{\left(x_{k-\mathrm{1}}-\mathit{\mu }\right)}^{\mathrm{2}}}{{\mathit{\sigma }}^{\mathrm{2}}}}+{\displaystyle \frac{{\left(x_k-\mathit{\mu }\right)}^{\mathrm{2}}}{{\mathit{\sigma }}^{\mathrm{2}}}}\\ & -{\displaystyle \frac{\mathrm{2}\mathrm{\varnothing }\left(x_{k-\mathrm{1}}-\mathit{\mu }\right)\left(x_k-\mathit{\mu }\right)}{{\mathit{\sigma }}^{\mathrm{2}}}}\left]\right).\end{array}\end{array}$$

Then the conditional probability distribution of X_t given $X_{t-\mathrm{1}}=c$ can be derived by the Bayes' theorem and written as (see Appendix A)

$$\begin{array}{}\text{(6)}& p\left(x_t\mathrm{|}{x}_{t-\mathrm{1}}=c\right)\sim N\left(\mathit{\mu }+\mathrm{\varnothing }\left(c-\mathit{\mu }\right),\left(\mathrm{1}-{\mathrm{\varnothing }}^{\mathrm{2}}\right){\mathit{\sigma }}^{\mathrm{2}}\right),\end{array}$$

where c is an arbitrary constant. The implication of such a formulation is that the resulting probability is also a function of c: if the statistical model parameters $(\mathit{\mu },{\mathit{\sigma }}^{\mathrm{2}},\mathrm{\varnothing })$ are fixed, a shorter distance of c from the mean, μ, will result in a relatively higher probability density than those are far away.

2.2 Constant value episode (CVE) probability

The estimation of the CVT probability consists of the following two steps.

Step 1. Deriving a joint probability density. For a series of (dependent) events, A_k with $\mathrm{1}\le k\le n$, the joint density of probability can be described through a product of multiple conditional probabilities as

$$\begin{array}{}\text{(7)}& \begin{array}{rl}p\left(A_n\cap \mathrm{\dots }\cap A_{\mathrm{1}}\right)& =p\left(A_{\mathrm{1}}\right){\prod }_{k=\mathrm{2}}^{n}p\left(A_k\mathrm{|}{\cap }_{j=\mathrm{1}}^{k-\mathrm{1}}{A}_j\right)\\ & =p\left(A_{\mathrm{1}}\right){\prod }_{k=\mathrm{2}}^{n}p\left(A_k\mathrm{|}{A}_{k-\mathrm{1}}\right).\end{array}\end{array}$$

The first equality yields from the chain rule of the joint distribution (Schum, 2001); the second equality is a special case of an AR(1) process.

Step 2. Imposing a distributional assumption to the joint probability distribution. From Eq. (6), the probability of consecutive values in a series with Gaussian probability density can be determined by

$$\begin{array}{}\text{(8)}& \begin{array}{rl}P\left({\text{CVE}}_{t=\mathrm{1},c\ne \mathrm{0}}\right)& =p(x_t=c\mathrm{|}{x}_{t-\mathrm{1}}=c)\\ & =\underset{c-\mathrm{res}/\mathrm{2}}{\overset{c+\mathrm{res}/\mathrm{2}}{\int }}{\displaystyle \frac{\mathrm{1}}{\mathit{\sigma }\sqrt{\mathrm{2}\mathit{\pi }\left(\mathrm{1}-{\mathrm{\varnothing }}^{\mathrm{2}}\right)}}}\\ & \exp \left(-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[{\displaystyle \frac{{\left(\left(c-\mathit{\mu }\right)-\mathrm{\varnothing }\left(c-\mathit{\mu }\right)\right)}^{\mathrm{2}}}{\left(\mathrm{1}-{\mathrm{\varnothing }}^{\mathrm{2}}\right){\mathit{\sigma }}^{\mathrm{2}}}}\right]\right).\end{array}\end{array}$$

The integral reflects the fact that digital data are recorded with finite numerical precision. Then, according to the property of an AR(1) process, the probability of a CVE with a length of t can be calculated through P(CVE₁) raising to the power of t−1 as

$$\begin{array}{}\text{(9)}& \begin{array}{rl}P\left({\text{CVE}}_{t,c\ne \mathrm{0}}\right)& =(\underset{c-\mathrm{res}/\mathrm{2}}{\overset{c+\mathrm{res}/\mathrm{2}}{\int }}{\displaystyle \frac{\mathrm{1}}{\mathit{\sigma }\sqrt{\mathrm{2}\mathit{\pi }\left(\mathrm{1}-{\mathrm{\varnothing }}^{\mathrm{2}}\right)}}}\\ & \exp \left(-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[{\displaystyle \frac{{\left(\left(c-\mathit{\mu }\right)-\mathrm{\varnothing }\left(c-\mathit{\mu }\right)\right)}^{\mathrm{2}}}{\left(\mathrm{1}-{\mathrm{\varnothing }}^{\mathrm{2}}\right){\mathit{\sigma }}^{\mathrm{2}}}}\right]\right){)}^{t-\mathrm{1}}.\end{array}\end{array}$$

Since this equation is designed for a constant event, so the marginal probability remains a constant for each CVE. To diminish the influence of CVE on μ, they were excluded first, then the μ, σ and ∅ were calculated.

For non-normal cases, the explicit parameterization of a non-independent joint distribution is difficult to derive due to mathematical challenges and often does not have a closed form. The non-parametric alternative is to use empirical distribution (Epanechnikov, 1969; Waterman and Whiteman, 1978) or kernel distribution (Hwang et al., 1994; Duong and Hazelton, 2005), but this approach is not desirable for database management at this stage because it is difficult to develop a unified framework that is adequate for all situations. Besides, the empirical distribution estimates a probability without taking into account auto-correlation, i.e., independent of the adjacent data points.

The AR(1) assumption can be relaxed by increasing the order of autocorrelation without too much complexity. For example, for an AR(2) process, one could specify the covariance matrix in Eq. (1) as

$$\begin{array}{}\text{(10)}& \mathrm{\Sigma }=\left|\begin{array}{ccc}{\mathit{\sigma }}^{\mathrm{2}}& {\mathit{\sigma }}^{\mathrm{2}}{\mathrm{\varnothing }}_{\mathrm{1}}& {\mathit{\sigma }}^{\mathrm{2}}{\mathrm{\varnothing }}_{\mathrm{2}}\\ {\mathit{\sigma }}^{\mathrm{2}}{\mathrm{\varnothing }}_{\mathrm{1}}& {\mathit{\sigma }}^{\mathrm{2}}& {\mathit{\sigma }}^{\mathrm{2}}{\mathrm{\varnothing }}_{\mathrm{1}}\\ {\mathit{\sigma }}^{\mathrm{2}}{\mathrm{\varnothing }}_{\mathrm{2}}& {\mathit{\sigma }}^{\mathrm{2}}{\mathrm{\varnothing }}_{\mathrm{1}}& {\mathit{\sigma }}^{\mathrm{2}}\end{array}\right|\end{array}$$

and modify Eq. (7) in step 1 as

$$\begin{array}{}\text{(11)}& \begin{array}{rl}& p\left(A_n\cap \mathrm{\dots }\cap A_{\mathrm{1}}\right)\\ & =p\left(A_{\mathrm{1}}\right)p\left(A_{\mathrm{2}}\mathrm{|}{A}_{\mathrm{1}}\right){\prod }_{k=\mathrm{3}}^{n}p\left(A_k\mathrm{|}{A}_{k-\mathrm{1}},A_{k-\mathrm{2}}\right),\end{array}\end{array}$$

then update the conditional probability parameterized by $(\mathit{\mu },{\mathit{\sigma }}^{\mathrm{2}},{\mathrm{\varnothing }}_{\mathrm{1}},{\mathrm{\varnothing }}_{\mathrm{2}})$ in step 2. The more general extension of the autoregressive model is out of the scope of this study and can be referred to in Box et al. (2015).

For the variables with extra incidences of zero such as nitrogen oxides (NO, NO₂) and ozone, the lower interval of the integration in Eq. (9) was changed from c− res to 0. Note that in reality “zero” values in measurements may actually be recorded as small positive or negative numbers. This detail is ignored in the following because there is no universally applicable correction available. Some data sets may require a linear or non-linear bias correction, while for other data sets a simple cutoff, e.g., set to zero if | value $|<$ threshold, may be more appropriate.

3 Model sensitivity test

The P in Eq. (9) is affected by the parameters μ, σ, ∅, c, t and res. A simulation study was developed to evaluate the sensitivity of P to each parameter. Several experiments were conducted by generating a synthetic data series to demonstrate the influence of each parameter. For each experiment, the CVT was performed over a range of possible values.

A set of first-order autoregressive, AR(1), time series with hourly time steps and a length of 240 values (10 d) was generated using Eq. (2) and a random noise generator. As a reference case (ref), we set μ=10, σ=4 and ∅=0.8. The numerical precision was defined as 0.01. Four sets of CVEs with the same length (t=3) were added to this time series. The distance of the CVE from the mean, i.e., c−μ, was given as 0, 1, 2 and 3σ (see Fig. 1). In this figure, four CVEs are illustrated with a color code, i.e., red, blue, cyan and black, which are shown with boxes. The P varies from $\mathrm{7.67}\times {\mathrm{10}}^{-\mathrm{6}}$ for the first CVE to $\mathrm{4.77}\times {\mathrm{10}}^{-\mathrm{7}}$ for the fourth (last) CVE. As stated in Sect. 2.1, the value of P decreases as c−μ increases. CVEs which are further away from the mean are less likely to occur in nature.

Figure 1 A synthetic AR(1) time series with Gaussian data distribution and four arbitrarily selected CVEs of length t=3, with μ=10, σ=4, ∅=0.8, and $c-\mathit{\mu }=\mathrm{0}$, 4, 8, and 12, respectively. The CVEs are shown using a color code, i.e., red, blue, cyan and black. The numerical precision (res) is chosen as 0.01.

To assess the effect of t on P, a set of values ranging from 2 to 10 were selected for the t. All other parameters were fixed as in the baseline time series. As expected from Eq. (9), the P decreases exponentially with t (Fig. B1a). Note that the slope of this exponential decrease depends on c−μ. The larger the c−μ, the larger would be the slope. That is in agreement with Fig. 1, where the P decreases as the CVE gets further from the mean. However, the probability of finding two consecutive data points with the same value is about 1:300, i.e., in a year-long time series such incidents are expected to occur naturally about once per year if the sampling resolution is daily and about 25 times if the sampling resolution is hourly.

To investigate the non-linear influence of σ on P in Eq. (9), a range of values, i.e., 0.1, 0.2, 0.3, 0.4, 0.5, 1, 2, 3, 4, 5, 10 and 20, were set as σ, while other parameters remained unchanged. In this scenario, the P changes from $\mathrm{1.22}\times {\mathrm{10}}^{-\mathrm{2}}$ for the smallest σ to $\mathrm{8.93}\times {\mathrm{10}}^{-\mathrm{8}}$ for the largest one (Fig. B1b). By using Eq. (9), it thus becomes possible to estimate likelihoods for naturally occurring CVEs for data sets with different variability, in contrast to classical approaches, which use a fixed variability threshold.

The most interesting parameter to consider in the CVT is the lag-1 auto-correlation (∅). A sensitivity experiment with several additional time series was performed to assess the sensitivity of P with respect to ∅ (Fig. B1c). In this figure, P ranges from $\mathrm{1.23}\times {\mathrm{10}}^{-\mathrm{10}}$ to $\mathrm{2.5}\times {\mathrm{10}}^{-\mathrm{3}}$. The larger the ∅ (i.e., stronger persistence), the larger would be the probability of naturally occurring CVEs. The estimated probability is very sensitive to ∅ as it approaches 1. At the limit value of 1, Eq. (9) is undefined. If ∅=0, the time series only consists of noise, so it is less probable to get any CVEs.

Another parameter influencing P is the data digital resolution (res) or precision, where the data have been recorded in a fixed numerical precision (number of decimals) or as integers with possible rounding to the nearest multiple of 5, 10, etc. This parameter is shown in Eq. (9), where the resulting probability is integrated over the range of values from c− res$/\mathrm{2}$ to c+ res$/\mathrm{2}$.

To investigate the sensitivity of the P to the res parameter, the baseline time series was resampled by using several resolutions, i.e., 0.0001, 0.0002, 0.0005, 0.001, 0.002, 0.005, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1, 2 and 5. As shown in Fig. B2a for the example of res =5, larger res leads to additional CVEs, and it becomes harder to distinguish valid episodes from erroneous incidents. But here, to isolate influence of res on P, first the data were truncated to a new resolution, then the CVEs were added to the data. The CVT results are shown in Fig. B2b, in which the P changes from $\mathrm{4.77}\times {\mathrm{10}}^{-\mathrm{11}}$ to $\mathrm{7.57}\times {\mathrm{10}}^{-\mathrm{1}}$. This shows that by increasing the res, the P increases, meaning that if the data are recorded in a coarse resolution, there is a higher chance to count those data as valid data.

An experiment with several scaling factors, i.e., fc =0.1, 0.2, 0.5, 1, 2, 5 and 10, was performed to check the robustness of the CVT to the different data transformations. In this experiment, the CVEs were added first; then the scaling, i.e., x(t)× fc, was applied; and the data were truncated to a new numerical resolution given by res × fc. Scaling changes other parameters such as μ or σ, except ∅, which remains invariant. Figure B1d shows the robustness of the CVT output (P) with scaling. It is important to note that Eq. (9) is robust to the other data transformation such as normalization and standardization (see Appendix C).

A combined sensitivity analysis was performed to illustrate the effect of the parameters σ, ∅ and res in Eq. (9), i.e., the conditional probability for two consecutive values was evaluated over a range of conditions (σ and ∅ from 0.01 to 0.99, and res of 0.01, 0.1 and 0.5), with $\mathit{\mu }-c=\mathrm{0}$. The results are shown in Fig. 2 and can be interpreted as an upper limit for P that two successive values are valid data because $\mathit{\mu }-c=\mathrm{0}$ represents the maximum of the Gaussian distribution in Eq. (9). Using the chain rule from Eq. (11), these results can easily be extrapolated to longer CVEs. As Fig. 2 shows, the probability of finding two valid consecutive data points with the same value decreases rather quickly with increasing standard deviation σ. The ∅ has limited influence up to values of around 0.7. Above this threshold, the likelihood of a two-value CVE increases drastically. A coarser numerical resolution makes it more likely to encounter constant values in reality. At res similar to σ, the length, t, of the CVE will have to be much larger than 2 to reliably classify it as erroneous. In practical applications, one would generally set a threshold for the acceptable probability first. The information provided in Fig. 2 can then help to identify typical parameters of the time series, where this threshold will be reached.

Figure 2 Conditional probabilities to find a measured value x_t given x_t−1 for three different numerical resolutions, i.e., (a) res =0.01, (b) res =0.1 and (c) res =0.5. In this figure, the σ and ∅ range from 0.01 to 0.99.

4 Results and discussion

Two data time series were retrieved from the Tropospheric Ozone Assessment Report (TOAR) database (Schultz et al., 2017) to illustrate the practical use of the CVT. This database holds in situ measured data time series for ground-level ozone in hourly time resolutions. We selected the time series of ozone mixing ratio at the Azusa station (34^∘8^′ N, 117^∘55^′ W) in California that has data from the 1980s, when the data were recorded with a resolution of 8 ppb. Besides this, the TOAR database contains data for meteorological variables at some stations. We selected one temperature time series at the Cape Grim station, Tasmania (40^∘68^′ S, 144^∘69^′ E). This station is located at an altitude of 94 m directly on the coast, and it is a Southern Hemisphere background site with an extensive record back into 1980. The station primarily measures air which has passed over the Southern Ocean for several days. So, temperature variations at this site are often of small amplitude. Data series of carbon monoxide at the Fresno station (36.78^∘ N, 119.77^∘ W) were obtained from the EPA AQS database. This data was reported with a precision of 1 ppm in 1980 and later in 2022 with a higher precision of 0.001 ppm. The precision changes might have arisen from the method detection limits (0.5 and 0.001), measurement methods (instrumental non-dispersive infrared and instrumental gas filter correlation Teledyne API 300 EU) or method types (non-FRM and FRM; Federal Reference Method) detailed in their data files.

4.1 Temperature

Temperature is one of the key variables relevant to air quality research. For example, temperature is often used as a primary predictor for smog-related air quality. For demonstration of the CVT in a real data situation, 10 d of a temperature time series were selected. The μ, σ and ∅ of the selected 10 d time series are 12.55, 1.59 and 0.94, respectively. The recorded numerical resolution of the data is 0.01. The time series along with the probability, P, of each value being a valid observation is shown in Fig. 3. Altogether, 18 CVEs are visible in Fig. 3, 15 of them with t=2, 2 with t=3 and 1 with t=4.

Figure 3 Temperature time series at the Cape Grim station (40^∘68^′ S, 144^∘69^′ E) from 10 to 20 January 1983. Black and blue lines show the temperature value (^∘C) and its associated probability, P, in Eq. (9), respectively. In this figure, the time is shown in UTC. The P is not affected by the unit conversion, i.e., degrees Celsius (^∘C) to kelvin (K). The data were retrieved from the TOAR database.

The CVEs occur at more or less regular times in the early morning, e.g., 04, 05 and nighttime hours, e.g., 10, 21, 22 and 23 (see Fig. 4). That can be because of the local meteorological phenomena at this site, where the temperature has little variance. Therefore, these CVEs are less likely to be erroneous data.

Figure 4 The number of CVEs occurring for the different hours in a day, i.e., h={0…23}, for the temperature time series shown in Fig. 3.

The probabilities estimated by the CVT are above 0.2 in most cases, which means that if the CVEs were to be flagged as erroneous data, one would err in one out of five cases and throw out the valid measurements. The CVE on 18 January yields the lowest probability (0.008), in line with the expectation of the human data analyst because it is a sparse CVE with four consecutive values (t=4). This example illustrates that it will generally be impossible to define a universal threshold for P but that instead it depends on the use case. For example, in a data quality control workflow at the originating institution, one may decide to rule out data with $P<{\mathrm{10}}^{-\mathrm{4}}$ but have a data curator cross-check the measurements with larger P. In contrast, when these data are integrated in a larger analysis consisting of many stations, one might apply the CVT to rule out data with $P<{\mathrm{10}}^{-\mathrm{3}}$ or even $P<{\mathrm{10}}^{-\mathrm{2}}$ to increase the statistical robustness of the analysis.

Other criteria for selecting a threshold for P could be climate regions. In the polar regions, the diurnal cycle of the temperature in summer could be quite high, but coastal sites in that area with a dense fog might have morning periods when the temperature is rather constant. The first shows a larger σ than the latter, so the P will be less in the polar than the coastal sites, assuming all other parameters are constant (as shown in Fig. B1b). One may adopt a smaller threshold for P in polar than in coastal sites. Or for the same climatological region with constant temperature values at night or in the day, when the diurnal cycle reaches maximum or minimum, the CVT would give CVEs a lower probability, as they are further from the mean (larger c−μ). So, the P of the CVEs at extrema can be less than the CVEs with the same t in this series.

4.2 Ozone

Ozone near the ground is an air pollutant that is detrimental to human health and vegetation growth. Ozone measurement techniques have evolved over time, and it can therefore be challenging to assess the data quality of a decade-long monitoring data set, such as that from the Azusa station in California, U.S. (34^∘8^′ N, 117^∘55^′ W), that contains a relatively long data record from 1980 to 2016.

Figure 5 shows a 10 d example from this measurement series for the year 1990, with the μ, σ and ∅ of 16.55, 17.32 and 0.79, respectively. During the early period, the data were reported in a low resolution, here an interval of 8 ppb. As a consequence, the time series contains many CVEs, and most of them are probably valid. In contrast, for the year 2012 when the data are recorded in a higher data resolution, i.e., 1 ppb, the number of the CVE is small (see Fig. D1). As mentioned in the introduction, urban ozone time series often show very low values (effectively zero), which are, however, recorded as small positive or negative values, here +2 ppb. Figure 5 shows the probabilities between $\mathrm{3.12}\times {\mathrm{10}}^{-\mathrm{10}}$ and 1 for these episodes, which have values of 2 ppb. There are also three CVEs, with large t (≥8) and very low ozone mixing ratios of 2 ppb, which are shown with red circles in Fig. 5. This illustrates the issue of zero-bounded data mentioned in the methodology. The CVT can recognize such cases, and the associated probabilities are $\mathrm{3.12}\times {\mathrm{10}}^{-\mathrm{10}}$, $\mathrm{2.22}\times {\mathrm{10}}^{-\mathrm{7}}$ and $\mathrm{2.48}\times {\mathrm{10}}^{-\mathrm{8}}$, for the CVE1, CVE2 and CVE3, respectively. That would prevent such (valid) values from being flagged or filtered as erroneous data, in contrast to the second part of the time series in Fig. 6 (for the year 2011), which exhibits sparse occurrence of episodes, i.e., 21 CVEs where 17, 2, 1 and 1 CVEs with the t=2, 4, 7 and 9, respectively. In most cases (17 episodes), the CVEs consist of only two consecutive values (t=2). The estimated probability for these cases is between $\mathrm{2.15}\times {\mathrm{10}}^{-\mathrm{2}}$ and $\mathrm{9.9}\times {\mathrm{10}}^{-\mathrm{2}}$ (Fig. 6). One episode during 18 November 2011 consists of nine constant values of 2 ppb. The estimated P for that incident is $\mathrm{4.6}\times {\mathrm{10}}^{-\mathrm{14}}$, and this episode would indeed raise the suspicions of trained data analysts because such a pattern in the data would require a rather special explanation (see Fig. D3).

Figure 5 Time series of the ozone mixing ratio at the Azusa station, California, from 10 to 20 November 1990 (black) and the CVT test results (blue). During this period, the data were recorded in intervals of 8 ppb, i.e., res =8, so that valid CVEs are frequent. In total, this time series contains 45 CVEs as 27, 6, 3, 3, 1, 1, 1 and 1 episode, with the t=2, 3, 4, 5, 6, 8, 9, and 11, respectively. The red circles (or ovals) highlight three examples of zero-ozone incidents (here 2 ppb) with a large length (t≥8) in this series. The cyan circles highlight the probability of the respective CVEs. The orange circle highlights a CVE with a length of 4 that contain a gap of missing data points.

Figure 6 As Fig. 5 but from 10 to 20 November 2011, when the data were recorded with a numerical resolution of 1 ppb, i.e., res =1. The red circle shows one example of missing data points in the data time series. The μ, σ and ∅ of the data in this figure are 19.9, 10.73 and 0.84, respectively.

Figure 5 also illustrates the problem with missing data values that was mentioned in the beginning of Sect. 2. On 18 November, there is a gap in the time series where the data point has been excluded, and the values to the left and right of this episode are identical. If these values were not treated correctly, they would be counted as a CVE episode with a length of 8 and probability of $\mathrm{2.58}\times {\mathrm{10}}^{-\mathrm{7}}$, which is shown with an orange circle in Fig. 5. Although such incidents could raise suspicions, they are not (and should not be) detected by the CVT. An independent test needs to be designed for such situations.

4.3 Carbon monoxide

Exposure to elevated carbon monoxide harms the human body, in particular those who suffer from heart diseases. This air pollutant also affects some greenhouse gases, e.g., carbon dioxide and ozone, which are linked to climate change and global warming. A 10 d example of the measured carbon monoxide at the Fresno station is shown in Fig. 7. Despite the high precision of the data for the year 2022 (res =0.001, see Fig. D4), data were recorded with a resolution of 1 ppm in 1980. These data contain fewer CVEs but with a larger t (19 CVEs with t=2…34) in comparison to the ozone series in Fig. 5. That could be associated with a longer lifetime of carbon monoxide than that of ozone. This reflects that most of the CVEs in the carbon monoxide series are valid. The CVT discerns this and estimates a larger P for this data, in which the smallest P is 0.001 for the CVEs, with t=14 and values of 0 ppm.

Figure 7 Time series of carbon monoxide at the Fresno station, California, from 1 to 11 January 1980 (black) and the CVT test results (blue). During this period, the data were recorded in intervals of 1 ppm, i.e., res =1, so that valid CVEs are frequent. In total, this time series contains 19 CVEs as 1, 1, 1, 1, 2, 2, 1, 2, 1, 1, 1, 3 and 2 episodes with the t=34, 27, 21, 18, 15, 14, 12, 11, 10, 5, 4, 3 and 2, respectively. The μ, σ and ∅ of the data in this figure are 0.79, 0.45 and 0.65, respectively.

5 Conclusions

Environmental time series are valuable and essential data sources for scientific assessment of air quality and climate change. One of the issues in these data is the occurrence of the constant value episodes (CVEs). These episodes are often considered to be indicative of sensors' malfunctions or other measurement errors and are excluded from the data via quality control (QC) procedures. However, these episodes can be due to the natural environmental phenomena, and they are indeed valid observations. Thus, distinguishing whether the CVEs are erroneous or valid data is accompanied by large uncertainty.

This study presented a theoretical concept and evaluation for a data-driven constant value test (CVT), which takes into account the typical evolution of environmental state variables such as air temperature, ozone mixing ratio or carbon monoxide as time series with serial dependence. Based on the calculus of a marginal, joint and conditional Gaussian probability density, one can estimate the probability of constant value episodes (CVEs) of length t to occur in reality and use this information to flag data as potentially erroneous. The threshold for such flagging needs to be selected by the data analyst. Together with the batch size for processing pieces of the time series (in our examples, the full length of the depicted data was used; for practical applications on longer time series, we recommend sample sizes in the order of 100), these are the only a priori parameters needed. Examples with synthetic and real data demonstrate that the CVT captures many aspects which a trained data analyst would consider in the QC of such time series. But as a data-driven approach, it will reveal data inconsistencies (here, CVEs due to measurement or data processing errors) in automated data processing workflows, and it may assist manual data quality control by making it possible to provide a fine-grained warning to the data analyst that something may be wrong with the measurements based on a probabilistic score.

The test first detects CVEs by testing for zero difference. Then, it evaluates the distribution parameters mean (μ), standard deviation (σ) and lag-1 auto-correlation (∅), as well as the numerical resolution of the data in user-defined portions (batches) of the time series. Given these parameters, the conditional probability for two consecutive identical values is computed and integrated over the interval given by the numerical resolution of the recorded data. Using the chain rule for the non-independent conditional probability, this probability can easily be scaled to arbitrary lengths of CVEs.

The novelty of this approach is its foundation in statistical theory and the concept of estimating the probability of a data sample to occur naturally. This distinguishes the method from classical approaches where more or less arbitrary thresholds need to be defined prior to testing. Such pre-defined thresholds can be dangerous if conditions change, for example, when the same thresholds are applied to data from different world regions, climatic zones or seasons. The method is robust against such changes, and its application requires little background knowledge about the specific data set under investigation. The method is therefore well suited for having robust and automated QC systems, for example, in smart sensor networks, where human intervention is not feasible.

Appendix A

The inference of conditional probability of bivariate normal distribution

$$\begin{array}{rl}& {\displaystyle \frac{f\left(x_{k-\mathrm{1}},x_k\right)}{f\left(x_{k-\mathrm{1}}\right)}}\\ & ={\displaystyle \frac{\frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }{\mathit{\sigma }}^{\mathrm{2}}\sqrt{\mathrm{1}-{\mathrm{\varnothing }}^{\mathrm{2}}}}\exp \left(-\frac{\mathrm{1}}{\mathrm{2}\left(\mathrm{1}-{\mathrm{\varnothing }}^{\mathrm{2}}\right)}\left[\frac{{\left(x_{k-\mathrm{1}}-\mathit{\mu }\right)}^{\mathrm{2}}}{{\mathit{\sigma }}^{\mathrm{2}}}+\frac{{\left(x_k-\mathit{\mu }\right)}^{\mathrm{2}}}{{\mathit{\sigma }}^{\mathrm{2}}}-\frac{\mathrm{2}\mathrm{\varnothing }\left(x_{k-\mathrm{1}}-\mathit{\mu }\right)\left(x_k-\mathit{\mu }\right)}{{\mathit{\sigma }}^{\mathrm{2}}}\right]\right)}{\frac{\mathrm{1}}{\mathit{\sigma }\sqrt{\mathrm{2}\mathit{\pi }}}\exp \left(-\frac{\mathrm{1}}{\mathrm{2}}\left[\frac{{\left(x_{k-\mathrm{1}}-\mathit{\mu }\right)}^{\mathrm{2}}}{{\mathit{\sigma }}^{\mathrm{2}}}\right]\right)}}\\ & ={\displaystyle \frac{\mathrm{1}}{\mathit{\sigma }\sqrt{\mathrm{2}\mathit{\pi }\left(\mathrm{1}-{\mathrm{\varnothing }}^{\mathrm{2}}\right)}}}\exp (-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}(\mathrm{1}-{\mathrm{\varnothing }}^{\mathrm{2}})}}[{\displaystyle \frac{{\left(x_{k-\mathrm{1}}-\mathit{\mu }\right)}^{\mathrm{2}}}{{\mathit{\sigma }}^{\mathrm{2}}}}+{\displaystyle \frac{{\left(x_k-\mathit{\mu }\right)}^{\mathrm{2}}}{{\mathit{\sigma }}^{\mathrm{2}}}}\\ & -{\displaystyle \frac{\mathrm{2}\mathrm{\varnothing }\left(x_{k-\mathrm{1}}-\mathit{\mu }\right)\left(x_k-\mathit{\mu }\right)}{{\mathit{\sigma }}^{\mathrm{2}}}}-{\displaystyle \frac{\left(\mathrm{1}-{\mathrm{\varnothing }}^{\mathrm{2}}\right){\left(x_{k-\mathrm{1}}-\mathit{\mu }\right)}^{\mathrm{2}}}{{\mathit{\sigma }}^{\mathrm{2}}}}\left]\right)\\ & ={\displaystyle \frac{\mathrm{1}}{\mathit{\sigma }\sqrt{\mathrm{2}\mathit{\pi }(\mathrm{1}-{\mathrm{\varnothing }}^{\mathrm{2}})}}}\exp (-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}\left(\mathrm{1}-{\mathrm{\varnothing }}^{\mathrm{2}}\right)}}[{\displaystyle \frac{{\mathrm{\varnothing }}^{\mathrm{2}}{\left(x_{k-\mathrm{1}}-\mathit{\mu }\right)}^{\mathrm{2}}}{{\mathit{\sigma }}^{\mathrm{2}}}}\\ & +{\displaystyle \frac{{\left(x_k-\mathit{\mu }\right)}^{\mathrm{2}}}{{\mathit{\sigma }}^{\mathrm{2}}}}-{\displaystyle \frac{\mathrm{2}\mathrm{\varnothing }\left(x_{k-\mathrm{1}}-\mathit{\mu }\right)\left(x_k-\mathit{\mu }\right)}{{\mathit{\sigma }}^{\mathrm{2}}}}\left]\right)\\ & \sim N(\mathit{\mu }+\mathrm{\varnothing }\left(c-\mathit{\mu }\right),(\mathrm{1}-{\mathrm{\varnothing }}^{\mathrm{2}}\left){\mathit{\sigma }}^{\mathrm{2}}\right),\text{given}{x}_{k-\mathrm{1}}=c.\end{array}$$

Appendix B

Figure B1 Sensitivity of P to the (a) CVEs length, i.e., t=2, 3, 4, 5, 6, 7, 8, 9 and 10. Other parameters are fixed as μ=10, σ=4, ∅=0.8 and $c-\mathit{\mu }=\mathrm{0}$, 4, 8 and 12. (b) Standard deviation, i.e., σ=0.1, 0.2, 0.3, 0.4, 0.5, 1, 2, 3, 4, 5, 10 and 20. Other parameters are fixed as μ=10, t=3, ∅=0.8 and $c-\mathit{\mu }=\mathrm{0}$, 4, 8 and 12. (c) Lag-1 autocorrelation, i.e., ∅=0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.91, 0.92, 0.93, 0.94, 0.5, 0.96, 0.97, 0.98 and 0.99. Other parameters are fixed as μ=10, σ=4, t=3 and $c-\mathit{\mu }=\mathrm{0}$, 4, 8 and 12. (d) Sensitivity of P to scaling factor, i.e., fc =0.1, 0.2, 0.5, 1, 2, 5 and 10. Other parameters are fixed as ∅=0.8 and t=3. The same color codes are applied as those in Fig. 1.

Figure B2 (a) The modified time series (res =5) where ref time series were resampled with rounding to the nearest of five. That includes more CVEs than the ref in Fig. 1. (b) Sensitivity of P to the digital numerical precision, i.e., res =0.0001, 0.0002, 0.0005, 0.001, 0.002, 0.005, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1, 2 and 5. Other parameters are fixed as μ=10, σ=4, ∅=0.8, t=3 and $c-\mathit{\mu }=\mathrm{0}$, 4, 8 and 12. The same color codes are applied as those in Fig. 1.

Appendix C

If the data are normalized, i.e., $(x-x_{min})/(x_{max}-x_{min})$

Figure C1 As Fig. 1 but the data time series are normalized, μ=0.5, σ=0.15, ∅=0.8 and res =0.004.

If the data are standardized, i.e., $(x-\mathit{\mu })/\mathit{\sigma }$

Figure C2 As Fig. 1 but the data time series are standardized, $\mathit{\mu }=-\mathrm{0.07}$, σ=0.94, ∅=0.8 and res =0.002.

Appendix D

Figure D1 Time series of the ozone mixing ratio at the Azusa station, California, from 10 to 20 November 2011. During this period, the data were recorded in intervals of 1 ppb, i.e., res =1. μ=19.9, σ=10.73 and ∅=0.84.

Figure D2 As Fig. 6, but the missing values are not treated. So, the orange circle shows two CVEs, which have been merged to one incident with a longer length (t=8).

Figure D3 Number of CVEs (∑t) of different length, i.e., $t=\left\{\mathrm{0},\mathrm{\dots },\mathrm{9}\right\}$, for the ozone time series of the year 2011 (shown in Fig. 6).

Figure D4 As Fig. 7 but from 1 to 11 January 2022, when the data were recorded with a numerical resolution of 0.001 ppm, i.e., res =0.001. This series shows three CVEs with the length of 2, i.e., t=2. The μ, σ and ∅ of the data in this figure are 0.62, 0.4 and 0.88, respectively.