Method for the Detection of Functional Outliers Applied to Quality Monitoring Samples in the Vicinity of El Musel Seaport in the Metropolitan Area of Gijón (Northern Spain)

Abstract

Air pollution affects human health and is one of the main problems in the world, including in coastal cities with industrial seaports. In this sense, the city of Gijón (northern Spain) stands out as one of the 20 Spanish cities with the worst air quality. The study aims to identify outliers in air quality observations near the El Musel seaport, resulting from the emissions of six pollutants over an eight-year period (2014–2021). It compares methods based on the functional data analysis (FDA) approach and vector methods to determine the optimal approach for detecting outliers and supporting air quality control. Our approach involves analyzing air pollutant observations as a set of curves rather than vectors. Therefore, in the FDA approach, curves are constructed to provide the best fit to isolated data points, resulting in a collection of continuous functions. These functions capture the behavior of the data in a continuous domain. Two FDA approach methodologies were used here: the functional bagplot and the high-density region (HDR) boxplot. Finally, outlier detection using the FDA approach was found to be more powerful than the vector methods and the functional bagplot method detected more outliers than the HDR boxplot.

1. Introduction

Air pollution is a global problem with significant health implications [1]. The growth in maritime transport due to globalization has made seaports crucial for economic and social development, accounting for over 80% of world trade [2,3]. Several ports in northwest Spain, including La Coruña, Ferrol, and Gijón (El Musel), have expanded to accommodate increased traffic [4,5]. Pollution issues have also been observed in other Spanish ports such as Santander and La Coruña, which have implemented various solutions to mitigate the problem [6].

Port activities, including ship maneuvers, fueling, and docking, contribute to substantial pollution, especially during non-optimal engine operation phases [7,8]. Ships, particularly those with diesel engines, emit significant amounts of particulate matter, both primary and secondary, which can have a substantial impact on nearby populations [9,10]. Industrial ports additionally contribute to air pollutant emissions through discharges, handling of bulk materials, road dust, and heavy vehicle traffic, all influenced by meteorological conditions, material types, and proximity to populations [11,12]. In some cases, traffic within the port area can have a more significant impact on air quality than bulk unloading [13].

All these activities have been identified as sources of air pollution that can negatively affect human health, leading to the exceeding of daily limit values for inhalable particles [14,15]. Numerous studies have examined the impact of air quality in cities with seaports and the transport of pollution inland [16,17]. Coastal areas, especially those downwind of port areas, experience specific meteorological patterns that affect the dispersion, transformation, accumulation, and deposition of air pollutants [18].

Several studies have evaluated the impact of loading and unloading operations on particulate matter emissions in ports, such as those conducted in the ports of Tarragona and Thessaloniki [12,19]. Other studies have examined PM₁₀ (particles smaller than 10 microns) and PM_2.5 (particles smaller than 2.5 microns) emissions resulting from port activities in Barcelona and Civitavecchia, as well as noise and pollution caused by ship and port activities in Ancora [20,21,22]. Pollution from ships in ports such as those in Shanghai and Istanbul has been analyzed, along with studies on pollutant gas compositions, PM₁₀ levels due to bulk handling in Alicante, and the relationship between maritime traffic and pollutant concentration in several ports [17,23].

It is well-known that there are different types of outliers, such as vector and functional outliers. To date, there have been a number of studies on the identification of vector outliers in fields as diverse as medicine [24] and renewable energy [25]. Numerous methods have been proposed to detect and investigate outliers in pollutant measurements recorded at air quality monitoring stations [26,27,28]. However, no single method has been universally accepted as the definitive approach.

Outlier detection based on vectorial methods, such as Z-score and Rosner test, involves analysis of the distribution of individual data points with respect to the mean and standard deviation of the data set. On the other hand, methods based on the FDA approach are designed to analyze functional data sets, where each observation is a curve or function. These outlier detection methods are Bagplot and HDR Boxplot.

In the area surrounding the port of Gijón, there have been recent and repeated episodes of pollution that have caused social alarm and concern about possible health effects. Municipal protocols to combat pollution have also been activated on several occasions. The objective of this study was to propose a model for the identification of outliers derived from pollutant emissions into the atmosphere in the vicinity of the industrial port of El Musel in Gijón (latitude 43°34′ N, longitude 5°41′ W) using different methodologies. These included vectorial methods and methods based on the FDA approach, allowing a comparative analysis between the two approaches. In this way, early warning indicators can be provided to help local authorities make timely decisions.

The structure of this paper is as follows: Section 2 introduces the experimental setup and variables involved in this study. Section 2.4 explains the mathematical methods used in this paper, including Box–Cox transformation, outlier detection using vector methods, smoothing technique, functional bagplot, and HDR boxplot. Section 3 describes the results obtained with the proposed models and provides a discussion by comparing both methodologies. Finally, Section 4 concludes this study by highlighting the main findings.

2. Materials and Methods

2.1. Gijón and the Port of El Musel

Gijón, with a population of 270,871 registered in 2021, is located on the coast of the Cantabrian Sea (Figure 1a), in the autonomous community of Principality of Asturias, in northern Spain. Gijón’s climate is characterized by short, mild summers with an average temperature of 18.7 degrees Celsius (°C) in August. Winters are long, wet, and relatively cold, with an average minimum temperature of 8.1 °C in January and February, which are the coldest months. It rains throughout the year in Gijón, with July being the driest month with an average of 72 millimeters (mm) and November being the wettest with an average of 159 mm. The total average annual rainfall is 1303 mm and the average annual temperature is 13.0 °C. Prevailing wind conditions are SW in winter and NE in summer. The strongest wind conditions are NW.

The construction of the Port of El Musel in Gijón was due to the impulse of the iron and steel industry and the exploitation of hard coal in the second quarter of the 19th century in Asturias. Commercial exploitation began in 1907 [29]. It has undergone numerous expansions throughout its history, the most recent of which took place between 2005 and 2011. Currently, El Musel is an artificial outer harbor separated from the urban area of the city of Gijón, whose main economic activity is the traffic of goods, connected to more than 200 ports worldwide. The port is composed of seven main piers. Figure 1b shows its location in relation to the city. The total length of quays and berths in the commercial basins are 6934 meters (m) and 1006 m in fishing basins. The main materials unloaded include coal, petroleum coke, clinker, iron ore, and other non-metallic materials. It is the largest bulk port in Spain, with 14,500 thousand tons (kt) in 2019 and 13,616 kt of solid bulk in 2020 [30,31]. This can cause a large amount of air pollution, which is further exacerbated by road traffic along the port facilities and the roads leading to and from the port.

2.2. Pollution

Air pollutants can be divided into two types. The first type is primary pollutants, which are emitted directly from sources into the atmosphere. The second type is secondary pollutants, which are formed from primary pollutants via chemical reactions and microphysical processes in the atmosphere [32].

The origin of pollutants can be attributed to human activities. In this case, they are called anthropogenic. Otherwise, they are referred to as naturally occurring pollutants.

Sulfur dioxide (SO₂) is produced, among other sources, by exhaust emissions from diesel engines, including those from ships, vehicles, and industrial plants [33]. Aerosols formed via reaction with water vapor are hazardous to health even at low concentrations and can cause irritation in the respiratory system [32]. They also cause damage to structures by contributing to corrosion problems through the formation of sulfuric acid and producing acid rain [34].

Carbon monoxide (CO) is a gas produced primarily by anthropogenic sources, such as incomplete combustion from road traffic [35].

Ozone (O₃) is a highly irritating pollutant to humans, is harmful to plants, and is commonly found in urban smog. It is of secondary origin, formed through the reaction of NO_x (nitrogen oxides) and VOCs (volatile organic compounds) with sunlight. It causes numerous respiratory and health problems, as reported in numerous studies [36,37].

PM10 particles are small enough to pass through the respiratory tract and lodge in the lungs. Short-term exposure can aggravate existing health issues, while long-term exposure may lead to lung disease and other ailments [38].

The most harmful nitrogen oxides as pollutants are NO₂ and NO (nitrogen dioxide and nitric oxide) [32]. Extensive research has shown the harmful effects of these pollutants on human health, affecting both adults [39] and children [40]. They can react with water vapor to form nitric acid, which, like sulfur oxides, causes metal corrosion issues [41].

2.3. Experimental Data Set

In this article, the values of the pollutants PM₁₀, NO, NO₂, SO₂, CO, and O₃ collected over six years, from January 2016 to December 2021, at the Principality of Asturias monitoring station ES1271A, located in Gijón, have been analyzed. This air quality monitoring station, called Argentina, is the closest to the port, as shown in Figure 1b.

2.4. Mathematical Methodologies Used in This Study

2.4.1. Box–Cox Transformation

In order to carry out the study of vector outliers, a normality condition was required, which was not present in the initial data recorded by the air quality monitoring station of this study. One of the objectives of variable transformation is to reduce skewness as well as other distributional properties that cause difficulties in statistical analysis. Therefore, it was necessary to find a simple transformation that could approximate the data of a variable to a normal distribution. Linear transformations do not change the shape measures of skewness and kurtosis [42]. However, nonlinear transformations can modify these measures and can be used to obtain variables with distributions that approximate normal distribution. Among these transformations, the most commonly used is the Box–Cox transformation, which belongs to the following family of transformations [43]:

$$x_j^{\left(\lambda \right)}=\left\{\begin{array}{c}\frac{x_j^{\lambda }-1}{\lambda }, if \lambda \ne 0\\ \log \left(x_j\right), if \lambda =0\end{array}\right.,$$

(1)

where $x_j$ represents the variable to be transformed and $\lambda$ is the transformation coefficient that maximizes the likelihood function of the variable’s profile. This makes it possible to perform a transformation that allows the data of a study to be treated as approximately normal in most cases [44].

Subsequently, after applying the Box–Cox transformation, normality was assessed using the Anderson–Darling test, a widely known test described in [45].

2.4.2. Outlier Detection Techniques Using Vector Methods

Z-Score

One methodology employed to identify atypical values (or outliers) is well-known as the Z-score, which assumes data normality [46]. This approach identifies outliers as standardized values that exceed a given threshold. The Z-score is calculated as follows:

$$z_i=\frac{x_i-x_r}{\sigma },$$

(2)

where $x_i$ represents the value of the observation at each point, $x_r$ is a reference value (typically the mean of the observations), and $\sigma$ is the standard deviation of all observations. The detailed development of this methodology can be found in [46].

The scoring criterion specified by the International Organization for Standardization/International Electrotechnical Commission (ISO/IEC 17043:2010) was adopted as the criterion to establish possible outliers [47]:

$$\begin{array}{c}\left|z_i\right|<2\to \mathrm{Satisfactory};\\ 2\le \left|z_i\right|\le 3\to \mathrm{Uncertain};\\ \left|z_i\right|>3\to \mathrm{Unsatisfactory}.\end{array}$$

(3)

Therefore, a Z-score greater than two was considered a possible outlier, while a score greater than three indicated a clear outlier. It is important to note that this criterion assumes a normal distribution of the contamination values, necessitating performing a prior transformation using the Box–Cox methodology.

A set of observations may have upper and lower outliers, which are, respectively, significantly higher or lower than the remaining values. In the case of air quality studies, higher concentrations of pollutants are of particular concern to human health, and therefore require monitoring. By applying the Z-score methodology, positive and negative values can be obtained. Following the criteria given in Equation (3) to identify possible outliers, we identified higher outliers ($z_i>3$) and lower outliers ($z_i<3$). However, for the purpose of this study, only positive values ($z_i>3$) were considered. Since our focus was on air quality, our interest was only in the upper outliers.

Rosner’s Test

Rosner’s test is a statistical tool used to detect outliers within a data set. Its basic principle is to compare the difference between a given value and the mean of the remaining values, which is then divided by the standard deviation of those values. A larger standardized difference indicates a higher probability that the value is an outlier.

The test is performed iteratively, meaning that the procedure is repeated for each potential outlier value, with the value with the largest standardized difference being eliminated at each iteration.

Furthermore, the specific method employed in Rosner’s test is known as the generalized ESD (Extreme Studentized Deviate) Test. This approach addresses certain limitations found in the Grubbs Test, also referred to as ESD or Pearson-Hartley. The ESD Test was developed to tackle two primary limitations: first, it was designed to test a single outlier and cannot find multiple outliers in a data set; second, it can suffer from a masking effect due to the presence of other outliers nearby [48].

Rosner’s test is designated to detect one or more outliers in a data set having an approximately normal distribution. It requires specifying an upper limit, denoted as $r$, for the number of outliers to be tested [49]. The generalized ESD performs $r$ separate tests, looking for from one to $r$ outliers. This approach avoids the masking effect and allows for the detection of multiple outliers [49].

Rosner’s test hypotheses are as follows [50]:

H₀: 

There are no outliers in the given data set, i.e., the data set is normally distributed.

H₁: 

There are up to a maximum of r outliers in the given data set.

The development of Rosner’s test, indicated in [49,51], is as follows:

Let $x_1,x_2,\dots ,x_n$ represent the $n$ observations of a given data set $X$. It is considered that $n-r$ of the observations of that data set come from a normal distribution, while the $r$ observations are potentially extreme and may or may not belong to a different distribution.

Let us now consider a new data set $X^{*}$ consisting of the elements $x_1^{*},x_2^{*},\dots ,x_{n-i}^{*}$, where $n-i$ are the number of observations after having removed the $i$ most extreme measurements from the previous set $X$; therefore, let $i=0,1,2,\dots ,r-1$.

We denote the standard deviation as $s^{\left(i\right)}$ and the mean as ${\stackrel{-}{x}}^{\left(i\right)}$ of the set $X^{*}$, as well as ${\stackrel{-}{x}}^{\left(0\right)}$ and $s^{\left(0\right)}$ to the mean and standard deviation of the set $X$ containing all observations.

We can define the most extreme observation in the data set, a specific value of i, as that which is furthest from the mean. That is:

$$x^{\left(i\right)}={max}_j\left|x_j^{*}-{\stackrel{-}{x}}^{\left(i\right)}\right|,j=1,2,\dots ,n-1.$$

(4)

Although, from the point of view of contaminant concentration analysis, the largest extreme observations are of interest, formally, from a mathematical point of view, the absolute value is taken to consider as an extreme observation the smallest or the largest of the elements of the data set.

Rosner’s outlier test [49,50] is supported by $r$ statistics, denoted as $R_1,R_2,\dots ,R_r;r\in \mathbb{N}∕r<n$, corresponding to the studentized extreme deviations obtained from subsets of decreasing size, i.e., $n,n-1,\dots ,n-r+1$ elements.

Therefore, we can define the statistic $R_{i+1}$ and the critical value ${\lambda }_{i+1}$ as follows:

$$R_{i+1}=\frac{\left|x^{\left(i\right)}-{\stackrel{-}{x}}^{\left(i\right)}\right|}{s^{\left(i\right)}},$$

(5)

$${\lambda }_{i+1}=\frac{t_{p,n-i-2}\left(n-i-1\right)}{\sqrt{\left(n-i-2+t_{p,n-i-2}\right)\left(n-i\right)}},$$

(6)

where $p=1-\frac{\alpha }{2\left(n-i\right)}$, $t_{p,\nu }$ is the $p^{th}$ quantile with $\nu$ degrees of freedom (dof) of the t-Student distribution, and $\alpha$ is the Type I error level.

The procedure to establish the number of outliers is as follows: first, compare the statistic $R_r$ with the critical value ${\lambda }_r$. If $R_r>{\lambda }_r$, then we can conclude that the most extreme $r$ value is an outlier. If $R_r\le {\lambda }_r$, proceed to compare $R_{r-1}$ with ${\lambda }_{r-1}$. If $R_{r-1}>{\lambda }_{r-1}$, then the $r-1$ most extreme values are considered outliers. Finally, we continue in the same way until we have detected up to a certain number $r$ of outliers or until the test does not give any value [50].

Rosner’s test is indicated for sets with $n>25$ observations [50].

2.4.3. Outlier Detection Methods Based on the FDA Approach

Functional data analysis (FDA) is a statistical framework designed for analyzing data in the form of smooth curves or functions. The basic philosophy of FDA is that, despite their discrete nature, observed data functions should be considered as single entities, rather than as a mere sequence of individual observations. A functional data set, denoted by $x$, consists of $n$ pairs $(t_j,Y_j)$, where $Y_j$ represents a recording or observation of the function $x\left(t_j\right)$ at argument value $t_j$. Since observation noise is typically present in data, the representation of raw data in a functional form often involves a smoothing process. The FDA approach is developed in detail in Ramsay and Silverman [52].

Outlier detection within the FDA approach focuses on identifying unusual observations or curves that deviate significantly from the expected behavior. These outliers provide valuable insight into abnormal or extreme patterns and are of great interest in various fields. The goals of outlier detection using methods based on the FDA approach include identifying outliers of magnitude and shape, characterizing their impact on the data, and understanding the underlying mechanisms that contribute to their occurrence.

The following subsections describe smoothing in more detail, as well as the two outlier detection methods based on the FDA approach used in this study: functional bagplot and functional high-density region boxplot.

Smoothing: A Method to Represent Curves Based on Points

In certain problem domains, measurements are derived from operations that have functional characteristics. In the present study, the data under investigation were derived from the average daily concentrations of atmospheric pollutants, aggregated on a monthly basis, as recorded by an air quality monitoring station. The focus was on the analysis of a set of air pollutants observed over a consecutive period of months. Initial exploratory analysis revealed a discernible pattern of variation with a sinusoidal profile. This observed variation could potentially be influenced by meteorological factors such as temperature, rainfall, humidity, wind direction, and wind speed, as well as various processes associated with nearby industrial and port activities. The FDA approach represents an appropriate method for investigating these types of issues.

In the FDA approach, curves are created to fit the isolated data points, transforming a vector sample into a functional sample. This results in a set of continuous functions instead of multiple point values.

Applying appropriate smoothing assumptions in data analysis allows for a greater amount of information to be obtained and contributes to flexible and robust results. The smoothing approach must accurately represent the desired characteristics.

Smoothing and curve fitting are different concepts. Curve fitting aims to find a curve that closely resembles the original data.

According to Ramsay and Silverman [52], monitoring a continuous process at isolated points leads to monitoring functional data. Let $F$ be a functional space and $x\left(t_j\right)$ be a set of $n_p$ observations, where $t_j$ represents each instant of time. These observations can be considered as the values taken by the function $x\left(t\right)\in \chi \subset F$ at specific points. For evaluation purposes, we assume $F=span\left\{{\varphi }_1,\dots ,{\varphi }_{n_b}\right\}$, where $\left\{{\varphi }_k\right\}$ is a set of basis functions, $k=\mathrm{1,2},\dots ,n_b$, and $n_b$ is the number of functions needed to form the basis of the functional space. Typically, spline or Fourier functions are used as basis functions. Although other types are possible [53], we considered a Fourier function:

$$x\left(t\right)=\sum _{k=1}^{n_b}{c}_k{\varphi }_k\left(t\right),$$

(7)

where ${\left\{c_k\right\}}_{k=1}^{n_b}$ are the function coefficients of $x\left(t\right)$ with respect to a collection of selected basis functions. The smoothing problem is based on finding the solution to the following regularization issue [53]:

$$\underset{x\in F}{\min }\sum _{j=1}^{n_p}{\left\{z_j-x\left(t_j\right)\right\}}^2+\lambda \mathsf{\Gamma }\left(x\right).$$

(8)

In this expression, $z_j=x\left(t_j\right)+{\epsilon }_j$ calculates the value of the function $x$ at each point $t_j$, where ${\epsilon }_j$ is a random perturbation with zero mean. $\mathsf{\Gamma }$ is an operator that penalizes problem difficulty and $\lambda$ is a regularization coefficient controlling the intensity.

A Fourier series can approximate periodic functions, both univariate and multivariate. Thus, the expansion in Expression (8) can be written as:

$$\underset{c}{\min }\left\{{\left(\mathit{z}-\mathit{\Phi }\mathit{c}\right)}^T\left(\mathit{z}-\mathit{\Phi }\mathit{c}\right)+{\lambda \mathit{c}}^T\mathit{R}\mathit{c}\right\},$$

(9)

where $\mathit{z}={\left(z_1.\dots ,z_{n_p}\right)}^T$ is the vector from measurements, ${\mathit{c}=\left(c_1.\dots ,c_{n_p}\right)}^T$ is the vector of the coefficients from the functional expansion, $\mathit{\Phi }$ is an $n_p\times n_p$ matrix with elements ${\mathit{\Phi }}_{jk}={\varphi }_k\left(t_j\right)$, and $\mathit{R}$ is an $n_b\times n_b$ matrix defined as:

$$R_{kl}={〈D^2{\varphi }_k,D^2{\varphi }_l〉}_{L_2\left(T\right)}={\int }_T^{}{D}^2{\varphi }_k\left(t\right)D^2{\varphi }_l\left(t\right)dt.$$

(10)

In the above expressions, $D^n{\varphi }_k\left(t\right)$ represents the $n$th-order differential operator applied to ${\varphi }_k$. The solution can be obtained as:

$$\mathit{c}={\left({\mathit{\Phi }}^{\mathit{t}}\mathit{\Phi }+\lambda \mathit{R}\right)}^{-1}{\mathit{\Phi }}^t\mathit{z}.$$

(11)

In the FDA approach, the data from the measurements not only consist of discrete points, but also belong to a smooth function. Typically, this function exhibits some level of smoothness and continuity, where consecutive monitored points $z_j$ and $z_{j+1}$ are close to each other. The analysis is often performed on a group of functional data sets sharing a common performance pattern. In this study, the periodicity was considered by incorporating monthly monitoring.

Functional Bagplot

The concept of functional bagplot extends from the bivariate bagplot developed in [54] that made use of the first two robust scores of the principal component method.

The bivariate bagplot is a variation of the univariate boxplot proposed by Tukey [55]. It represents the semi-space position depth of an element $\theta \in {\mathbb{R}}^2$ with respect to a bivariate collection of points $Z=\left\{{\mathit{z}}_1,{\mathit{z}}_2,\dots ,{\mathit{z}}_n\right\}$. The localization depth of the semi-space $ldepth\left(\theta ,Z\right)$, $\theta \in {\mathbb{R}}^2$, $Z=\left\{{\mathit{z}}_1,{\mathit{z}}_2,\dots ,{\mathit{z}}_n\right\}$ is defined as the minimum number of ${\mathit{z}}_i$ included in any closed semi-plane with a boundary line passing through $\theta$. The depth zone $D_k$ is the collection of all points $\theta \in {\mathbb{R}}^2$, such as $ldepth\left(\theta ,Z\right)>k;k>0$. It is necessary to consider that the depth zones are all convex polygons satisfying $D_{k+1}\subset D_k$. The median depth of $Z$ is defined as the $\theta$ with the largest $ldepth\left(\theta ,Z\right)$.

Both the boxplot and the bivariate bagplot have the following components: a central element indicating the median depth, an inner case area, and an outer boundary area. Outliers are individual elements outside these areas. In a bagplot, the element with the highest depth (median depth) is located at the center, and it is surrounded by a case formed by 50% of the elements with the highest depth. This concept is analogous to the interquartile range in a traditional boxplot. By extending the case by a factor, the boundary is determined and any element outside this edge is considered an outlier. The bagplot not only detects outliers, but also provides information about correlation, position, as well as dispersion, skewness, tails, and a confidence zone for the depth median.

Functional bagplots extend the bagplot concept to a collection of curves $\left\{y_i\right\},$ $i=\mathrm{1,2},\dots ,n$ by utilizing the first two scores obtained through the robust principal components method of $y_i$ [56]:

$$y_i=\mu \left(x\right)+\sum _{k=1}^2{z}_{i,k}{\phi }_k\left(x\right),$$

(12)

where $\mu \left(x\right)$ represents the mean curve, $\left\{{\varphi }_k\left(x\right)\right\}$ are the principal components, and $z_i=\left(z_{i,1},z_{i,2}\right)$ are the first and second principal component scores. Analyzing the first two principal components enhances the detection of outliers compared with the original functional space [57]. These components typically capture the most important modes of variation [58]. The functional bagplot is particularly useful when the number of variables exceeds the amount of data monitored.

The functional bagplot presents the bivariate bagplot relative to the first and second scores. It includes the median curve, regions with 95% confidence for the median, the 50% case area, and 95% boundary area classified based on depth. Curves outside the 95% convex area are considered functional outliers. Outliers are positioned relative to the boundary in the bagplot.

A detailed development of the bagplot can be found in the study conducted by Ruts and Rousseeuw [59]. The key steps of this procedure are as follows:

• Calculate the number of elements included in the depth area $D_k$, denoted as $\#D_k$. Find the index $k$ for which $\#D_k\le \left(\frac{n}{2}\right)\le \#D_{k+1}$.
• Perform a linear interpolation, denoted as $f\left(x\right)$, based on the center of mass of the internal area labeled $T^{*}$ to obtain the case, which is a convex polygon in all cases.
• Determine the edge by increasing it with respect to $T^{*}$ by an index of 3. The choice of this quantity is supported by previous findings [54].
• Mark elements beyond the edge as outliers.

Functional High-Density Region Boxplot

In some studies [56], high-density regions (HDR) are considered a powerful tool to summarize prediction distributions. HDR exhibits versatility by combining modularity and asymmetry in the predictive density. It has the property of being the smallest feasible zone of probability coverage in a sample space at all times. In this article, we focus on a two-dimensional scenario which implies dealing with the minimum surface area.

The boxplot methodology [55] is very useful for studying univariate samples. Hyndman proposed a boxplot methodology based on HDR, which outlines probability distributions and can be extended to n-dimensional multivariate densities [60]. The HDR functional boxplot [56] is a modification of the functional boxplot that arranges the scores ${\mathit{z}}_i=\left(z_{i,1},z_{i,2}\right)$ using a bivariate kernel density approximation.

Let $f$ be a density function and $\left\{{\mathit{z}}_i;i=1,\dots n\right\}$ be a random sample of two variables from $f$. The kernel density approximation for Parzen–Rosenblatt in a bivariate case is as follows [61]:

$$\widehat{f}\left(\mathit{\omega };a,b\right)=\frac{1}{nab}\sum _{k=1}^{n}k\left(\frac{{\omega }_1-z_{i,1}}{a}\right)k\left(\frac{{\omega }_2-z_{i,2}}{b}\right),$$

(13)

where $\mathit{\omega }={\left({\omega }_1,{\omega }_2\right)}^{\prime };k\left(·\right)$ is a one-dimensional kernel function with symmetry, such that it is verified that $\int k\left(u\right)du=1$, and $\left(a,b\right)$, where $a>0$, $b>0$, control the bivariate bandwidth. The HDR boxplot can be defined as:

$$R_{\alpha }=\left\{\mathit{z}:\widehat{f}\left(z;a,b\right)\ge f_{\alpha }\right\},$$

(14)

where $f_{\alpha }$ ensures that $Pr\left(Z\in R_{\alpha }\right)\ge 1-\alpha$. The density estimated within this zone is greater than that outside the zone, with a coverage probability of $\left(1-\alpha \right)$. Notably, the HDR boxplot captures the smallest area in the sample space among all possible areas with a probability coverage of $\left(1-\alpha \right)$, and it includes the mode.

The functional HDR boxplot displays the functional outliers, as well as the mode, determined as ${sup}_z\widehat{f}\left(\mathit{z};a,b\right)$, and the inner and outer areas. The inner area corresponds to the region within the curves of the inner elements of the 50% bivariate HDR boxplot, akin to the inner quartile range of a boxplot. It provides information about the central 50% dispersion of the curves. The outer zone lies outside the curves of the inner elements of the 95% bivariate HDR boxplot and the outliers are the elements not contained within the 95% HDR boxplot.

To address complications arising from numerical integration of $f\left(x\right)$, Hyndman proposed a Monte Carlo technique [60]. This involves considering a multivariate random variable $X$ with density $f\left(x\right)$ and transforming it into $Y=f\left(X\right)$. A new function $f_{\alpha }$ is obtained by transforming $X$ using the density function, such that $Pr\left(f\left(X\right)\ge f_{\alpha }\right)=1-\alpha$ holds. Thus, $f_{\alpha }$ represents the quantile $\alpha$ of $Y$, and it can be estimated through a sampling process. More details can be found in Reference [62].

3. Results and Discussion

Air quality in the vicinity of the port of El Musel and the western area of Gijón can have a significant impact on people’s health and pose a risk factor for various diseases. In this study, we focus on analyzing pollutants based on data collected at the Argentina station, specially targeting pollutants with legislated values for health protection. These pollutants include:

• SO₂—Coal and oil often contain sulfur compounds and their combustion produces sulfur dioxide. Further oxidation of SO₂, usually in the presence of a catalyst such as NO₂, forms H₂SO₄ (sulfuric acid), resulting in acid rain. This is one of the reasons for concern about the environmental impact of the use of these fuels for power generation [32,63]. The limit values set for the protection of human health by European Community Directive 2008/50 [64] are 350 μg/m³ over an average period of 1 h and 125 μg/m³ over a 24 h mean period.
• NO₂—This gas is considered one of the most important atmospheric pollutants. It is present in the formation of acid rain as HNO₃ (nitric acid) due to reactions with water vapor droplets in the atmosphere [32]. The limit values set in the European legislation [64] are 240 μg/m³ over an average period of 1 h and 40 μg/m³ for the calendar year.
• CO—Carbon monoxide is a product of the incomplete combustion of fuels such as coal and natural gas. Automobile exhaust is also a source of carbon monoxide [32]. The maximum allowable concentration in air quality legislation [64] is 10 mg per day, based on a daily 8 h mean.
• PM₁₀—Particulate matter in the air represents a complex mixture of organic and inorganic substances. Particulate matter is found in urban areas from thermal power plants, industrial processes, vehicle traffic, domestic coal combustion, and industrial processes and activities [20]. The limit values for PM₁₀ collected in [64] are 40 μg/m³ per calendar year and 50 μg/m³ per day. This daily limit value may not be exceeded more than 35 times per calendar year.
• O₃—Tropospheric ozone of artificial origin is a secondary pollutant, as it is generated from other pollutants in the air due to a series of photochemical reactions of volatile organic compounds and nitrogen oxides [32]. The target value set by [64] legislation for the protection of human health is 120 μg/m³ as the daily maximum of the eight-hourly moving averages. This value shall not be exceeded for more than 25 days each calendar year considering a 3 year average.

Our main objective was focused on the analysis of daily ${\mathrm{P}\mathrm{M}}_{10}$ measurements, but daily measurements of ${\mathrm{S}\mathrm{O}}_2, \mathrm{N}\mathrm{O}, {\mathrm{N}\mathrm{O}}_2, \mathrm{C}\mathrm{O}, \text{and} {\mathrm{O}}_3$ pollutant gases from the same station were also studied. Since our objective was to extract and analyze statistical information from the collected data, it could happen that the identification of an outlier did not correspond to the consideration of exceeding the limits established by the legislation, since the latter values are fixed and the extreme values come from mathematical analysis and are therefore inherent in the data themselves.

For the pollutants studied, our sample ${\left\{x_{ij}^{\left(k\right)}\right\}}_{j=1}^n$ corresponded to $n=96$ months, spanning from January 2014 to December 2021. In this sample, $x_{ij}$ represents the concentration of each pollutant $k$ measured on day $i$ of month $j$, with $i=1,2,\dots ,30$ and $j=1,2,\dots ,n$, where $k=\left\{{\mathrm{S}\mathrm{O}}_2, \mathrm{N}\mathrm{O}, {\mathrm{N}\mathrm{O}}_2, \mathrm{C}\mathrm{O}, {\mathrm{P}\mathrm{M}}_{10}, {\mathrm{O}}_3\right\}$.

The results obtained for each of the methods used in the analysis of air pollution at the Argentina station, which was selected to characterize the emissions of the Gijón port area due to its proximity to the El Musel port area in Gijón, are presented below.

Table 1. Summary of the outliers detected using vector and functional procedures.

Method/Pollutant	SO2	NO	NO2	CO	PM10	O3
Vector			February 2015		March 2014	September 2018
					October 2014
					February 2020
HDR Boxplot	January 2015	November 2014	October 2014	October 2014	Mach 2014	May 2018
	November 2016	January 2015	February 2015	May 2017	October 2014	September 2018
	December 2016	December 2015	November 2018	November 2017	March 2015	May 2019
	November 2017	December 2016	January 2019	December 2018	February 2020
	December 2018	November 2018	March 2020	February 2019	December 2021
Bagplot	March 2014	October 2014	October 2014	October 2014	Mach 2014	May 2018
	January 2015	November 2014	February 2015	February 2015	October 2014	July 2018
	April 2015	February 2015	November 2018	December 2015	February 2015	September 2018
	February 2016	October 2015	January 2019	December 2016	March 2015	May 2019
	November 2016	November 2015	Mach 2020	April 2017	February 2020	June 2020
	December 2016	December 2015		May 2017	Mach 2021
	February 2017	January 2016		November 2017	December 2021
	April 2017	December 2017		December 2018
	November 2017	February 2018		February 2019
	November 2018	November 2018		October 2021
		December 2018
		February 2019
		January 2020
		December 2020

presents a summary of the outliers identified by the vector analysis and the two methods based on the FDA approach used in this study. It is important to note that the lower outliers were excluded from the analysis results, as they were not considered representative in the context of air pollution and air quality assessment.

Table 2. Month with maximum values as average concentration *.

Pollutant/Year	2014	2015	2016	2017	2018	2019	2020	2021
SO2 (µg/m3)	February (10)	January (17)	December (12)	February (12)	December (13)	February (12)	January (8)	December (8)
NO (µg/m3)	November (18)	December (30)	December (30)	January (20)	December (25)	December (16)	January (19)	December (17)
NO2 (µg/m3)	November (34)	January (38)	December (40)	January (37)	November (37)	February (38)	January (31)	December/January (27)
CO (mg/m3)	October (0.64)	December (0.68)	December (0.75)	May (0.86)	December (0.83)	February (0.93)	October (0.65)	January (0.50)
PM10 (µg/m3)	October (47)	December (44)	December (36)	October (35)	December (40)	February (42)	February (46)	December (49)
O3 (µg/m3)	July/May (46)	July (46)	April (47)	April (54)	September (65)	April (58)	April (53)	June/April (53)

* The average concentration is given in brackets. Source: adapted from [65].

presents the months and their respective monthly mean values, indicating the months in each year that recorded the highest monthly mean value.

3.1. Vector Analysis Results

During the data exploration phase, it was observed that the data set did not follow a normal distribution. An Anderson–Darling test was performed, yielding a p-value of less than 0.005 for all pollutants. Therefore, the normality of the data was not accepted at the 5% significance level.

A Box–Cox transformation was applied to the data to obtain a normal distribution. This transformation was appropriate because the pollutant values were always positive and there were no observations of pollutants equal to zero during the study period.

Table 3. Optimal values of λ for Box–Cox Transformation for SO₂, NO, NO₂, CO, PM₁₀, and O₃ pollutants.

Pollutant	Lower Confidence Limit	Upper Confidence Limit	Optimal λ
SO2	−0.42	−0.29	−0.35
NO	−0.11	−0.02	−0.07
NO2	0.26	0.40	0.33
CO	−0.44	−0.31	−0.38
PM10	0.06	0.23	0.15
O3	0.78	0.94	0.86

shows the optimal $\lambda$ value used for the transformation along the upper and lower limits of the 95% confidence interval for each pollutant.

The generalized ESD (extreme deviation studied) test was conducted as described above. A value of r = 1 for all pollutants, and for the pollutant PM₁₀, the test was repeated with r = 2 after having detected one outlier with r = 1. The results obtained are presented in

Table 4. Results of generalized ESD test.

Pollutant	i	Mean (i)	Standard Deviation (i)	Value	Observation (Day)	R (i + 1)	λ (i + 1)	Outlier
SO2	0	0.57981	0.11723	1.00000	2673	3.58428	4.28987	FALSE
NO	0	0.89639	0.04345	0.76782	665	2.95901	4.28987	FALSE
NO2	0	2.80590	0.45107	4.28506	1365	3.27923	4.28987	FALSE
CO	0	1.58098	0.35121	0.69273	1796	2.52911	4.28987	FALSE
PM10	0	1.62088	0.09149	2.02692	2611	4.43804	4.28987	TRUE
	1	1.62074	0.09119	1.26132	446	3.94134	4.28979	FALSE
O3	0	21.92891	7.14035	46.24131	1709	3.40493	4.28987	FALSE

. The Rosner test identified only one outlier for PM₁₀ in this study.

3.2. Results of the FDA Approach

By applying the described smoothing technique to the initial sample ${\left\{x_{ij}^{\left(k\right)}\right\}}_{j=1}^n$, a new sample $\left\{x_j\left(t\right)\right\}$ was generated for each analyzed pollutant $k$. In this new sample $x_j$ represents a new basis function derived from Fourier basis functions. In this analysis, a set of 100-element Fourier basis functions was specifically utilized. The correlation coefficient between the discrete values of the initial sample and the values of the functions generated from the new sample was 0.99, indicating a strong 99% correlation between the functional sample and the discrete values.

The FDA approach was performed after the samples were constructed using the functional smoothing method. Specifically, the two functional methods mentioned in the previous section were applied to a data set obtained from the Argentina air quality monitoring station.

Figure 2 and Figure 3 show the functional bagplot and the HDR boxplot, respectively, for the data collected at the Argentina station.

Both the functional bagplot and the HDR boxplot methods provide statistical means to identify outliers. Outliers are identified as curves that lie outside the light gray shaded areas that represent the bands corresponding to 95% confidence.

Table 1 shows the outliers obtained using the functional bagplot and HDR boxplot methods. Each outlier indicates the month in which a particular pollutant had a concentration outside the 95% confidence bands during the 8-year study period. The lower outliers were excluded from Table 1 as they were not considered representative for the purpose of the pollutant concentration analysis.

The highest values were recorded in the autumn and winter months for SO₂, NO, NO₂, CO, and PM₁₀ air pollutants. This is in agreement with observations made in several studies [66,67,68]. However, some exceptions can be seen, as in 2017 when the month with the highest monthly average concentration for the pollutant CO was May, with a value of 0.86 mg/m³, which corresponded to the highest value in the series after February 2019, with a value of 0.93 mg/m³.

By observing the graph of the daily concentration averages for the pollutant CO corresponding to the month of May for each of the 8 years studied (2014 to 2021) in Figure 4, it can be observed that the month of May presented an atypical behavior with values much higher than the rest of the annual series for the same month.

On the other hand, the pollutant O₃ presented the highest monthly mean values in spring and summer, as observed in some studies as [69], which correspond to the periods of highest solar activity. However, an exception can be considered in 2018, in which the month with the highest average concentration was September.

Regarding climatology, September 2018, as shown in Figure 5, was the third month with the most hours of sunshine with 183 h, compared with 251 in August and 214 in May. However, September 2018 was the month with the lowest precipitation month in 2018, with only 50.3 mm, following an atypical December with only 45.1 mm. The high sunshine hour activity combined with low precipitation in September 2018 could explain why the annual maximum of monthly averages was delayed until September in that year.

With vector analysis using the Z-score method, no outliers were detected for SO₂, NO, and CO pollutants. Only one outlier was detected for NO₂ corresponding to February 2015; three for PM₁₀ in March 2014, October 2014, and February 2020; and one outlier for O₃ in September 2018. Note that lower outliers were not considered as they were not representative in this air pollution study. The number of lower outliers detected using this method was PM₁₀ (five) and NO₂ (two). The Rosner test only detected one outlier for PM₁₀ in March 2021.

The number of outliers detected using the HDR boxplot and bagplot methods of the FDA approach was considerable, especially in the case of the bagplot. Outliers were detected for all contaminants using the FDA approach methods. Table 1 shows the outliers obtained by the vector and methods based on the FDA approach.

Based on the results obtained, it was observed that the outliers detected using vector analysis were also identified with the methods based on the FDA approach. This observation suggests that non-normality of the recorded measurements’ data set may be the reason for their presence. Classical methods are highly robust when dealing with data from normal distributions; however, their efficiency decreases when confronted with non-normal distributions.

Vector methods consider point observations to study the data. However, because air pollution from gases does not depend on single measurements, but on a series of measurements that may exceed the limits set by legislation within a specific time interval or, in this case, the statistically established limits derived from the data set, the FDA approach treats these observations as functional data rather than individual points. As a result, our approach involved treating air pollutant gas measurements as a collection of curves rather than as vectors, which can better capture variability. This perspective allowed us to view the original data set as a time-dependent function along a continuous domain rather than a set of discrete values at different points in time.

Methods based on the FDA approach are also more sensitive to detecting outliers than vector methods, especially when the data have complex or non-linear patterns. The FDA approach allows individual curves to be fitted to the data, which can be useful for detecting specific outliers in each curve.

However, they have some limitations, such as the fact that methods based on the FDA approach may require a higher level of knowledge and expertise compared with classical vector methods, coupled with higher computational complexity, which may increase the time and resources needed to perform the analysis. In the situation of data sets that follow a normal distribution, classical methods are generally quite robust and simpler. However, such classical methods lose efficiency on data sets that do not follow a normal distribution.

The proposed methods can be sensitive to the choice of the functional representation of the data and to the dimension reduction techniques. The choice of appropriate basis functions or dimension reduction methods can affect the quality of the functional representation and the subsequent detection of outliers.

Bias may also be introduced by the presence of incomplete or missing data. The methods assume complete and reliable observations, but in the real world, data may be missing for a variety of reasons, such as sensor failure or maintenance issues. Missing data can introduce bias and affect the accuracy of outlier detection.

As mentioned above, one advantage of methods based on the FDA approach over the identification of outliers by classical, i.e., vector methods, is the absence of need to assume normality. Therefore, there is no requirement to test for normality using methods such as Anderson–Darling or to perform a Box–Cox transformation.

In this study, the HDR boxplot method performed worse in detecting outliers of both methods based on the FDA approach. All outliers identified by the HDR boxplot method were detected by the bagplot method, except for two values for the pollutant NO, corresponding to months 13 (January 2015) and 36 (December 2016).

Figure 6 displays the graphs depicting the functions obtained after smoothing the pollutant emission data for SO₂, NO, NO₂, CO, PM₁₀, and O₃. Each pollutant had 96 functions representing each month of the 8-year study period conducted at the air quality monitoring station in Gijón.

With regard to maritime traffic, as indicated in the internal report of the Gijon Port Authority [70], the Port of El Musel had a significant decrease in the movement of goods through its facilities in 2020. It was a year marked by the health crisis and, in the case of the pier of Gijón, by the strong impact of decarbonization, as well as by the productive adjustment of steel production in the Asturian region. A slight recovery can be seen in 2021.

Table 5. Port maritime traffic in Gijón Port (2010–2021).

Year	2010	2011	2012	2013	2014	2015	2016	2017	2018	2019	2020	2021
Number of ships	1050	1065	1125	1251	1315	1183	1156	1298	1229	1187	1127	1189
Gross Tonnage												16,569
(thousands of tons)	14,254	17,735	19,091	20,504	21,192	19,097	18,311	19,607	18,226	18,005	16,468	9246
Dry Bulk												2021
(thousands of tons)	10,602	10,256	12,024	1669	12,382	15,355	12,036	14,768	12,717	9790	7010	1189
Number of ships	1050	1065	1125	1251	1315	1183	1156	1298	1229	1187	1127	16,569
Gross Tonnage												9246

Source: adapted from [31,71].

shows the evolution of maritime traffic in terms of number of vessels, gross tonnage, and dry bulk from 2010 to 2021.

Moldanová et al. [72] presented a PM emission factor from diesel-powered ships of 5.3 g, with a bimodal particle size distribution, one of the maxima being in the diameter of 7 microns and the other at 0.5 microns. Another study found that the total PM emission was 36 tons/year at Ambarli Port [73] and other researchers estimated 57.4 tons of PM per year for 7520 ships in the Candali Gulf [74], both located in Turkey.

According to Mueller et al. [75], factors such as road traffic, industrial and residential emissions, as well as ships’ emissions can affect air quality in port cities. The Argentina station in Gijón, due to its location, can receive not only emissions from the port of El Musel, but also emissions from the industrial areas to the west of Gijón and emissions from traffic in the area, as well as residential emissions. Another important factor is the daily truck traffic generated by the Port of El Musel. According to the consulting firm Vectio Traffic & Transport Planning (https://www.vectio.es) in a report prepared for the City Council of Gijón, the average number of trucks passing through one of the main access roads to the port was 1052 per day.

For the entire period 2004 to 2021, the Argentina station recorded exceedances of the limit values only in the case of the pollutant PM₁₀.

Table 6. Annual average of PM10 at the air quality station Argentina and number of exceedance days during the period 2004–2021.

Year	04	05	06	07	08	09	10	11	12	13	14	15	16	17	18	19	20	21
Exceedance days	76	93	84	70	25	45	41	63	44	41	36 *	21	6	14	16	11	26	38 *
Annual average (µg/m3)	48	43	42	39	31	36	33	37	34	35	33	28	25	29	31	28	31	38

(*) Due to discounts for Saharan intrusions, there were no instances of legislative non-compliance regarding exceedance days in 2014 and 2021. Exceeded values are highlighted in red. It is important to note that the current legislative values were revised in 2004 and have remained unchanged since then.

shows the number of exceedance days in that period and the annual average concentration of PM₁₀. The concentration of the others air pollutants were below the legal limits. Between the years 2009 and 2013, more than 35 exceedances per year of the limit value for the daily average concentration of PM₁₀, which is set at 50 µg/m³, were reached. In addition, there is the possibility of receiving pollutant emissions from several sources, as indicated above. These reasons have led to the implementation of an action plan to improve air quality.

For these reasons, the Government of the Principality of Asturias and the Gijón City Council commissioned a study from the Carlos III Institute of Health, which serves as the National Reference Laboratory for Air Quality. This study [76] recommended a new station located further west to better characterize industrial emissions.

Since the implementation of the actions in the 2014 Air Quality Improvement Plan by the Government of the Principality of Asturias and the City Council of Gijón, there has been a continuous decrease in the number of daily exceedances of PM₁₀ at the Argentina station (see Table 6), a trend that was broken in 2017 and 2018 with 14 and 16 daily exceedances, and subsequently in 2020 and 2021 with 26 and 38 exceedances, respectively. These data should be taken with caution, since at the end of this study period, the economic activities of shipping and industry were paralyzed or slowed down due to COVID-19, but there was an increase in the annual average PM₁₀ values, as well as in the number of exceedance days.

In 2021, the last year with complete data for this study, the Argentina station was in compliance with the number of allowable days exceeding the daily limit for PM₁₀, after applying discounts for Saharan intrusions. Saharan intrusions are a source of dust from the Sahara Desert [77].

Table 7. Average annual value and standard deviation * of each pollutant at the Argentina air quality station during the period 2004–2021.

Year	SO2 (µg/m3)	NO (µg/m3)	NO2 (µg/m3)	CO (mg/m3)	PM10 (µg/m3)	O3 (µg/m3)
2014	7.06 (1.99)	9.92 (5.31)	25.68 (6.95)	0.35 (0.12)	33.01 (7.49)	38.60 (8.72)
2015	7.23 (3.71)	12.21 (8.63)	27.29 (8.00)	0.35 (0.16)	28.38 (6.81)	35.75 (10.15)
2016	7.13 (2.10)	9.50 (7.79)	21.68 (8.48)	0.39 (0.16)	25.30 (4.64)	36.28 (9.55)
2017	8.11 (2.75)	9.27 (5.43)	26.25 (7.39)	0.50 (0.15)	28.85 (3.77)	38.92 (8.70)
2018	6.51 (3.30)	9.45 (7.20)	25.13 (7.53)	0.43 (0.17)	30.45 (4.81)	48.66 (12.31)
2019	5.49 (2.26)	10.54 (3.64)	22.75 (8.17)	0.48 (0.15)	27.73 (5.48)	44.25 (8.93)
2020	4.01 (1.76)	7.39 (5.20)	18.21 (7.15)	0.38 (0.15)	31.13 (9.26)	39.72 (10.18)
2021	4.37 (1.62)	8.47 (4.16)	18.53 (6.76)	0.27 (0.12)	38.06 (6.41)	40.56 (9.60)

(*) The value without brackets indicates the annual mean value and the value in brackets indicates the standard deviation.

shows the annual mean and standard deviation for each of the pollutants for the 8-year study period (2014–2021) at the Argentina air quality station.

As for the other air pollutants, all of them showed a decrease in 2020 compared with the values obtained in 2019. In 2021, the last year of the study, the annual average concentration of SO₂ increased slightly from 4 to 5 µg/m³. The same happened for NO, which increased from 7 to 9 µg/m³; with NO₂, from 18 to 19 µg/m³; and O₃ remained constant at 40 µg/m³ from 2020 to 2021. In contrast, CO continued the downward trend from 0.38 mg/m³ in 2020 to 0.27 mg/m³ in 2021. In 2021, there were no exceedances of the hourly limit value for any of the studied air pollutants.

The FDA approach for outlier detection has been used in other works, although with different approaches, either in the methods based on the FDA approach or in the field of application, such as those Díaz Muñiz et al. [78] applied to water quality control, in the emissions of polluting gases affecting urban areas [53]. Sguera et al. studied NOx levels in an industrial area to determine the variation in pollution between working days and weekends or holidays [28]. Chiou et al. [79] used an FDA approach to detect outliers in traffic flow data. In another study, the HDR boxplot and bagplot methods were used for the detection of outliers in data with hydrological applications [80]. In a recent study, new methods based on an FDA approach were proposed for identifying outliers in functional data sets. They were applied to climatological data series in Spain [81].

The calculations were performed on an AMD Ryzen 7 5800H processor, 3.20 GHz and 32 GB of RAM. The computational cost, using PM₁₀ as an example, was 4.01 s for the HDR and 0.08 s for the bagplot, compared with 0.02 s for the Z-score and less than 0.01 for the Rosner test.

4. Conclusions

In this study, we analyzed a set of data from measurements recorded by the air quality station of the Government of the Principality of Asturias, located closest to the Port of El Musel (Gijón), in order to detect outliers. This analysis was carried out because of the social alarm caused by the pollution episodes, especially PM₁₀, detected in this area, as an aid to decision making by local authorities, and for a possible means of early warning.

We used a vector method included in the European standards, such as the Z-score, a test to detect outliers, such as the generalized ESD test and the Rosner test, and more recent techniques to treat the discrete values of the observations as a series of continuous functions as a function of time.

Outlier detection using methods based on the FDA approach, such as HDR boxplot and bagplot, was more powerful than vector methods, such as Z-score and the generalized ESD test.

On the one hand, the analysis of the observations collected between 2014 and 2021 has shown the presence of the following upper outliers of air concentration for the vector method: NO₂ (one outlier), PM₁₀ (three), and O₃ (one). No outliers were detected for the pollutants SO₂, NO, and CO. The generalized ESD test detected a single outlier for PM₁₀ and none for the other pollutants.

The functional bagplot method detected a greater number of outliers than the HDR boxplot, with the following number of outliers being higher for the HDR boxplot method: SO₂ (five), NO (five), NO₂ (five), CO (five), PM₁₀ (five), and O₃ (three). The bagplot method identified SO₂ (10), NO (13), NO₂ (5), CO (10), PM₁₀ (7), and O₃ (5).

According to the results obtained in the pollutant analysis, the methods based on the FDA approach allowed us to retrieve more information than the vector methods. The latter compare averages, but do not take into account temporal variations.

Methods based on the FDA approach, in addition to having detected a greater number of outliers, have the advantage of not requiring that the data come from a normal distribution and, therefore, no hypothesis was made about the distribution followed by the recorded observations. This avoided the need to perform transformations on the data. In addition, this methodology considers the temporal correlation structure.

As for the origin of the pollution at the Argentina station, it could have several sources. In addition to the port of El Musel, it may have an industrial and traffic origin. For this reason, a study commissioned by the City Council of Gijón and the Government of the Principality of Asturias determined the need for a new station in a more western location. Therefore, new studies incorporating observations from this new station are needed to determine the origin.

Future research work, in addition to considering this new station, should focus on the application of this FDA approach methodology to statistical study and introduce sweep rules to evaluate trends in the functional samples.

In general, some of the factors that cause pollution to increase in an area are economic and population growth. In times of economic recession, pollution tends to decrease. It is possible to apply this air pollution outlier detection methodology to other air quality stations and in other geographic areas with success. It is also possible to apply this technique to a wide range of other types of problems, not necessarily related to air pollution, but in other areas such as water quality in water bodies, etc., to detect outliers. However, the specific characteristics of the site must always be taken in account.

This functional outlier detection methodology can be used as an early warning tool that allows authorities to take preventive decisions and effective measures to preserve the health of citizens.