A new approach for determining optimal placement of PM_2.5 air quality sensors: case study for the contiguous United States

We use a dataset consisting of daily PM_2.5 concentrations for the contiguous United States for the period January 2000 to December 2016. The dataset was produced through a data fusion method involving ensemble machine learning that combines surface monitoring measurements, satellite aerosol optical depth, land-use data, and chemical transport model results, among other variables (Di et al 2019, 2021). The high spatial and temporal resolution of this PM_2.5 dataset reveals spatial variability across the United States, including a stark east–west gradient, with the eastern United States experiencing relatively higher PM_2.5 concentration than the west. Hotspots in the dataset reveal the wide diversity of PM_2.5 sources, such as emissions from power plants in the Midwest, organic aerosol produced from biogenic volatile organic compounds in the southeast, and wildfires and associated smoke in the west. We coarsen the dataset's 1 km $\times$ 1 km resolution to 10 km $\times \,$ 10 km to yield a spatial grid of 871 $\times \,$ 413 cells. Missing values and erroneous observations are removed from the dataset by applying a mask to the spatial grid.

To isolate the underlying multiscale variability of PM_2.5, we use a DMD modeling approach. The DMD method provides a spatiotemporal decomposition of data into a set of dynamic modes derived from snapshots or measurements of a given system in time. DMD is similar to principal components analysis (PCA) as they are both dimensionality reduction algorithms that take advantage of underlying low-rank features found in large datasets. The resulting modes of these methods reveal the dominant spatial patterns that account for variation in the dataset. However, unlike PCA, which does not model temporal dynamics of time-series data explicitly, DMD computes a set of modes with an exponential temporal dependence whose imaginary components model an oscillating time frequency characterizing the modal dynamics. Importantly, DMD can accommodate heterogenous spatial and temporal sampling. For noisy datasets, DMD can be ensembled in order to produce uncertainty metrics (Sashidhar and Kutz 2021). Many variants of DMD exist (Tissot et al 2014, Noack et al 2015, Askham and Kutz 2017, Li et al 2017, Sashidhar and Kutz 2021), but here we follow Manohar et al (2018), as described in the supplement section 1 (available online at stacks.iop.org/ERL/17/034034/mmedia). We apply the DMD method to the PM_2.5 dataset to construct a library of modes that reveals the dominant spatiotemporal patterns in the dataset.

Recently, a multiresolution extension to the DMD algorithm has been proposed (Kutz et al 2016, Manohar et al 2019), which allows for the incorporation of significant transient phenomenon into the DMD library. Much like a multiresolution analysis using wavelets (Kutz 2013), the mrDMD algorithm recursively removes low frequency, or slowly varying, signals to allow separation of background (i.e. average PM_2.5 concentrations) from foreground (i.e. high pollution episodes) in the PM_2.5 dataset on varying timescales specified by the user. The mrDMD algorithm is more precise in capturing spatiotemporal variability than methods based on singular value decomposition such as DMD and PCA (Manohar et al 2019).

The number of snapshots M from our dataset are chosen so that the DMD modes provide a full rank approximation of the dynamics observed. Thus, M should be representative of all high- and low-frequency content of the dataset. The DMD algorithm is then applied recursively on the dataset. The slowest m1 modes are removed in the initial pass, and remaining snapshots are divided into M/2 segments and are fed into the DMD algorithm again. The slowest m2 modes are removed, and this process is continued until a termination condition is reached. The mrDMD separates the DMD approximate solution in the first pass as follows:

$$\begin{align}{X_{{\text{mrDMD}}}}\left( t \right) &= \sum\limits_{k{\,}={\,}1}^M {{b_k}} \left( 0 \right)\psi _k^{\left( 1 \right)}\left( \xi \right)\exp \left( {{\omega _k}t} \right)\nonumber\\& = \begin{array}{*{20}{c}} {\mathop \sum \limits_{k{\,}={\,}1}^{m1} {b_k}\left( 0 \right)\psi _k^{\left( 1 \right)}\left( \xi \right)\exp \left( {{\omega _k}t} \right)} \\ {\left( {{\text{slow modes}}} \right)} \end{array}\nonumber\\&\quad + \begin{array}{*{20}{c}} {\mathop \sum \limits_{k{\,}={\,}{m_1} + 1}^M {b_k}\left( 0 \right)\psi _k^{\left( 1 \right)}\left( \xi \right)\exp \left( {{\omega _k}t} \right)} \\ {\left( {{\text{fast modes}}} \right)} \end{array}\end{align} \tag{ 1 }$$

where ${b_k}\left( 0 \right)$ is the initial amplitude, $\xi$ is the spatial coordinates of each mode, $\psi _k^{\left( 1 \right)}$ represents the DMD modes (eigenvectors) obtained from the complete snapshot matrix, and $\exp \left( {{\omega _k}t} \right)$ are the corresponding eigenvalues. After the first level of decomposition, the fast modes are again decomposed into two matrices:

$$\begin{equation}{X_{M/2}} = X_{M/2}^{\left( 1 \right)} + X_{M/2}^{\left( 2 \right)}.\end{equation} \tag{ 2 }$$

At this level of decomposition, the first matrix $X_{M/2}^{\left( 1 \right)}$ represents slower modes and $X_{M/2}^{\left( 2 \right)}$ represents faster modes, which undergo decomposition again. This process is repeated recursively by removing slower modes in each iteration. At each level of decomposition, the separation of the foreground mode from the background is deemed significant if its eigenvalues exceed a user-set tolerance. Here, we set this tolerance to a standard value of 1 $\times$ 10⁻².

A formal expansion of the mrDMD theory and modeling approach may be found elsewhere (Kutz et al 2016, Manohar et al 2019).

The mrDMD algorithm operates by decreasing the time domain by a factor of two at each successive decomposition level. We apply mrDMD to training windows starting at 11.4 years (M= 4096 d), followed by 13 decomposition levels so that the shortest frequency is daily. This approach thus yields a long-term mode characterizing the average PM_2.5 concentrations over 11.4 years, with potential identification of transient pollution events spanning timescales from 5.7 years (2048 d) to 1 d. To analyze all PM_2.5 data from 2000 to 2016 employing these decomposition levels, we divide the dataset into two time-windows: January 2000 to March 2011 (first 4096 d) and September 2005 to December 2016 (last 4096 d).

We use the matrix libraries containing all DMD and mrDMD modes as tailored basis sets ${\psi _r} \in \,{\mathbb{R}^{n \times r}}$ to optimize for sensor placement. We identify the optimal sensor locations by employing QR pivoting to our DMD and mrDMD basis sets (Heck et al 1998, Manohar et al 2018). QR pivoting is a 'greedy' selection algorithm that is computationally efficient for finding near-optimal sensor locations. Greedy approaches are often favored over other optimization techniques as the true optimal solution often involves a combinatorially intractable optimization.

QR column pivoting identifies rows in the modal library ${\psi _r}$ with the highest two-norm, which corresponds to locations with the largest PM_2.5 modal frequencies and therefore greatest variability. The reduced matrix QR factorization with column pivoting decomposes a matrix $A \in \,{\mathbb{R}^{m \times n}}$ into a unitary matrix $Q$, an upper-triangular matrix $R$, and a column permutation matrix ${C^{\text{T}}}$ such that $A{C^{\text{T}}} = { }QR$. Thus, the QR factorization with column pivoting yields $r$ point sensors (pivots) that best sample the $r$ tailored basis modes ${\psi _r}$:

$$\begin{equation}\psi _r^{\text{T}}{C^{\text{T}}} = QR.\end{equation} \tag{ 3 }$$

That is, each QR pivot identifies those spatial locations in the modal library that exhibit the most variability and where sensor placement captures significant pollution episodes above background concentrations.

A formal expansion of the sparse sensor placement approach may be found elsewhere (Kutz et al 2016, Manohar et al 2018, 2019, Brunton and Kutz 2019).

In the supplement, we reconstruct the full field of PM_2.5 concentrations across the contiguous United States using measurements from only those sensor locations diagnosed as optimal in the step above. These reconstructions serve as validation of the DMD and mrDMD libraries and the ability of these algorithms to select sensor locations with high information content.

The mrDMD modal library reveals that most significant air pollution episodes occur on timescales between 1 week and 1 month, with many of these events arising from wildfires. We diagnose these wildfire modes by visual inspection based on regions with known fire activity and minimal industrial PM_2.5 sources. Figure 1 presents the resulting mrDMD modal maps in the time-frequency domain. Figure 1(a) corresponds to the first 11.4 years of the time period (January 2000 to March 2011) and figure 1(b) corresponds to the last 11.4 years (September 2005 to December 2016). The background mode in both training windows corresponds to the mean PM_2.5 concentration distribution for that window. Both background modes display a sharp east–west gradient with most of the eastern United States, except the Appalachian Mountains and some remote areas of Maine, having relatively higher PM_2.5 concentration than the west. This pattern agrees with the observed annual mean spatial distribution of PM_2.5 (Tai et al 2010, Di et al 2019). The background mode for the second half of the training period exhibits smaller modal amplitudes in hotspot regions than the first half. This trend is consistent with Di et al who reported that PM_2.5 concentrations decreased noticeably after 2008 due to a combination of the economic recession and emission controls on coal-fired power plants. The 2000–2008 annual average PM_2.5 concentration across the contiguous United States was 7.8 μg m⁻³; the 2009–2016 annual average was 6.2 μg m⁻³.

Zoom In Zoom Out Reset image size

The mrDMD method also identifies relatively short-term air pollution episodes on timescales of less than a year. Most episodes occur on timescales between 1 week and 1 month. However, we find an especially long-lived mode over Los Angeles and the San Joaquin Valley, California, lasting for 256 d from September 2007 to May 2008. This mode may correspond to intense winter PM_2.5 due to temperature inversions trapping emissions in the southern California basin. From 2007 to 2008, the San Joaquin Valley reported significantly higher PM_2.5 concentrations than are typical, with 66 d above the 24 h 35 μg m⁻³ standard in 2008 (US EPA 2021). The 2007 Murphy Complex Fire and other wildfires in Idaho are captured as a significant mode occurring on a timescale of 64 d from June to August. Although the lifetime of PM_2.5 in the atmosphere is ∼1 week and is sensitive to meteorology such as rainout, no pollution modes occur on timescales less than a week. We attribute the absence of modes on this timescale to the eigenvalue tolerance used in the mrDMD algorithm. If we tighten the tolerance, then more significant modes would be found at daily timescales, but such a tolerance would result in higher-frequency (i.e. noisy) and less-interpretable modal maps, making it more difficult to discern prolonged pollution events. Overall, we detect 72 significant pollution modes during the first time window and 77 during the second time window. We find 27 wildfire modes (38% of all modes) in the first window and 46 wildfire modes (60% of all modes) in the second window.

Based on figure 1, we conclude that the mrDMD framework can distinguish high pollution episodes from the mean PM_2.5 concentrations through time. The mrDMD framework provides a background mode that agrees with Di et al (2019), but also contains modal snapshots of events that stand out from this background.

By applying QR pivots to our mrDMD modes, we can identify a set of 1369 optimal sensor locations across the contiguous United States, a number of roughly the same magnitude as the number of EPA sensors. Our method assumes that air pollution spatiotemporal patterns and concentration variability from 2000 to 2016 are representative of future conditions. We argue that this assumption is reasonable, as discussed below. Figure 2 shows that the mrDMD sensor distribution is more clustered in California and across the west, along the Gulf Coast, and in the Industrial Midwest and Northeast, relative to that of the EPA. The mrDMD sensors thus reflect those regions characterized by higher PM_2.5 variability, with few sensors in the central United States. These results suggest that the west is a driver of PM_2.5 variability in the United States due to periodic wildfires, distributed urban sources isolated by rugged terrain, and reduced meteorological mixing. For example, the mrDMD algorithm yields a cluster of sensors in Los Angeles and Salt Lake City as these urban areas are basins surrounded by mountains where the meteorological conditions favor the production and trapping of PM_2.5 pollution. Additionally, we present a sparse distribution of sensors that can capture the overall spatiotemporal variability of PM_2.5 in the supplement.

Zoom In Zoom Out Reset image size

The mrDMD method identifies significantly more sensors in the San Joaquin Valley and in northern California, and the Pacific Northwest than the current EPA sensor network (figure 3). In particular, the San Joaquin Valley dominates the modal signatures of the mrDMD library (figure 1) and is characteristic of high PM_2.5 variability due to large-scale agricultural activities, smoke from wildfires, and frequent temperature inversions (Wolyn and McKee 1989, Wei et al 2013, Rémy et al 2017). The greater density of mrDMD sensors in northern California and the Pacific Northwest (especially in Idaho and the eastern parts of Oregon and Washington) can be attributed to their sensitivity to smoke PM_2.5 from wildfires. Indeed, we find that wildfires constitute 60% (46 out of 77) of significant pollution modes in the mrDMD library from September 2005 to December 2016. Also, mrDMD identifies locations along the Gulf Coast (i.e. the eastern shoreline of Texas) that are not covered by the EPA network (figure S1).

Zoom In Zoom Out Reset image size

Our results suggest that the EPA monitoring network does not optimally capture the spatial patterns and variability of PM_2.5 pollution in the contiguous United States today. We posit that the current locations of EPA sensors were informed by the legacy of air pollution from the 1980s to 1990s when the eastern United States had a higher regional PM_2.5 background with large variability in PM_2.5 concentrations arising from interactions between meteorology and emissions from coal-fired power plants, diesel emissions, and industrial activities (Bachmann 2007). However, with the passage of the Clean Air Act (and subsequent amendments), stricter air quality regulations have reduced these emissions (Bachmann 2007). Urban air pollution levels in turn declined nationally while at the same time the frequency and magnitude of wildfires in the west have increased significantly (Abatzoglou and Williams 2016). Our mrDMD sensor locations underscores this shift to the west as 53% of sensors (726 out of 1369) are found west of the 100th meridian (Seager et al 2017) compared to 32% (694 out of 2156) for the current network of EPA monitors. We emphasize that the current distribution EPA monitors is needed to keep track of adherence to the NAAQS, especially in densely populated regions. We do not advocate for the removal of existing PM_2.5 sensors, but rather for the deployment of more sensors in the west (especially in the San Joaquin Valley and Pacific Northwest) to capture high-pollution variability.