Robust Principal Component Analysis for Modal Decomposition of Corrupt Fluid Flows

TL;DR

This paper demonstrates that robust principal component analysis (RPCA) effectively cleans corrupted fluid flow data, preserving coherent structures and improving modal decomposition accuracy across simulations and experiments.

Contribution

It introduces the application of RPCA for robust modal analysis in fluid flows, handling corrupted data to enhance flow pattern identification.

Findings

RPCA accurately recovers low-rank flow structures.

RPCA effectively identifies and fills missing or incorrect data points.

RPCA improves the robustness of dynamic mode decomposition in noisy conditions.

Abstract

Modal analysis techniques are used to identify patterns and develop reduced-order models in a variety of fluid applications. However, experimentally acquired flow fields may be corrupted with incorrect and missing entries, which may degrade modal decomposition. Here we use robust principal component analysis (RPCA) to improve the quality of flow field data by leveraging global coherent structures to identify and replace spurious data points. RPCA is a robust variant of principal component analysis (PCA), also known as proper orthogonal decomposition (POD) in fluids, that decomposes a data matrix into the sum of a low-rank matrix containing coherent structures and a sparse matrix of outliers and corrupt entries. We apply RPCA filtering to a range of fluid simulations and experiments of varying complexities and assess the accuracy of low-rank structure recovery. First, we analyze direct numerical simulations of flow past a circular cylinder at Reynolds number 100 with artificial outliers, alongside similar PIV measurements at Reynolds number 413. Next, we apply RPCA filtering to a turbulent channel flow simulation from the Johns Hopkins Turbulence database, demonstrating that dominant coherent structures are preserved in the low-rank matrix. Finally, we investigate PIV measurements behind a two-bladed cross-flow turbine that exhibits both broadband and coherent phenomena. In all cases, we find that RPCA filtering extracts dominant coherent structures and identifies and fills in incorrect or missing measurements. The performance is particularly striking when flow fields are analyzed using dynamic mode decomposition, which is sensitive to noise and outliers.

Equations

$$\mathbf{X}$$

$$\mathbf{x}\approx \sum _{k=1}^r{\mathbf{u}}_k{\alpha }_k,$$

$${\mathbf{x}}_{k+1}\approx \mathbf{A}{\mathbf{x}}_k.$$

$$\mathbf{X}$$

$${\mathbf{X}}^{\prime }\approx \mathbf{AX}$$

$$\mathbf{A}=\underset{\mathbf{A}}{\mathrm{argmin}}\Vert {\mathbf{X}}^{\prime }-\mathbf{AX}{\Vert }_F={\mathbf{X}}^{\prime }{\mathbf{X}}^{†}\approx {\mathbf{X}}^{\prime }{\mathbf{V}}_r{\mathbf{\Sigma }}_r^{-1}{\mathbf{U}}_r^T.$$

$$\tilde{\mathbf{A}}={\mathbf{U}}_r^T\mathbf{A}{\mathbf{U}}_r={\mathbf{U}}_r^T{\mathbf{X}}^{\prime }{\mathbf{V}}_r{\mathbf{\Sigma }}_r^{-1}.$$

$$\tilde{\mathbf{A}}\mathbf{W}=\mathbf{W}\mathbf{\Gamma }.$$

$$\mathbf{\Phi }={\mathbf{X}}^{\prime }{\mathbf{V}}_r{\mathbf{\Sigma }}_r^{-1}\mathbf{W}.$$

$${\mathbf{x}}_k=\sum _{j=1}^r{\mathbf{\varphi }}_j{\gamma }^{k-1}{b}_j=\mathbf{\Phi }{\mathbf{\Gamma }}^{k-1}\mathbf{b},$$

$$\mathbf{b}={\mathbf{\Phi }}^{†}{\mathbf{x}}_1.$$

$$\mathbf{x}(t)=\sum _{j=1}^r{\mathbf{\varphi }}_j{e}^{{\omega }_{j}t}{b}_j=\mathbf{\Phi }\exp (\mathbf{\Omega }t)\mathbf{b},$$

$$\pm i\omega -ϵC{\omega }^2+\mathcal{O}({ϵ}^2),$$

$$\mathbf{X}=\mathbf{L}+\mathbf{S}.$$

$$\underset{\mathbf{L},\mathbf{S}}{\min }\text{rank}(\mathbf{L})+\Vert \mathbf{S}{\Vert }_0\text{subject to}\mathbf{L}+\mathbf{S}=\mathbf{X}.$$

$$\underset{\mathbf{L},\mathbf{S}}{\min }\Vert \mathbf{L}{\Vert }_{\ast }+{\lambda }_0\Vert \mathbf{S}{\Vert }_1\text{subject to}\mathbf{L}+\mathbf{S}=\mathbf{X},$$

$$\mathcal{L}(\mathbf{L},\mathbf{S},\mathbf{Y})=\Vert \mathbf{L}{\Vert }_{\ast }+{\lambda }_0\Vert \mathbf{S}{\Vert }_1+⟨\mathbf{Y},\mathbf{X}-\mathbf{L}-\mathbf{S}⟩+\frac{\upsilon }{2}\Vert \mathbf{X}-\mathbf{L}-\mathbf{S}{\Vert }_F^2.$$

$${\mathbf{Y}}_{k+1}={\mathbf{Y}}_k+\upsilon (\mathbf{X}-{\mathbf{L}}_k-{\mathbf{S}}_k),$$

$$\mathbf{L}$$

Full text

††∗ Corresponding author (ischerl@uw.edu).

Matlab code: github.com/ischerl/RPCA-PIV

Videos: tinyurl.com/RPCA-PIV

1 Introduction

The ability to understand, model, and control fluid flows is foundational to advancing technologies in transportation, energy, health, and defense. The challenges these fields pose are not easily solved by first principles analysis without intense simplification. Instead, we rely on data from simulations and experiments (taira2017aiaa, ; Duraisamy2018arfm, ; Brunton2020arfm, ). The fidelity of both approaches have improved dramatically, generating vast and increasing volumes of data Pollard2016book . However, despite this apparently high-dimensional data, fluid dynamics are often characterized by the evolution of a few dominant coherent structures that are energetically or dynamically important Lumley:1970 ; aubry1988dynamics ; berkooz1993proper ; HLBR_turb ; Kutz2016book . Thus, even with increasing ambient measurement dimension, the intrinsic dimension of the flow may remain relatively low. Modal decomposition techniques are designed to extract these meaningful patterns from high-dimensional fluids data taira2017aiaa ; taira2019aiaa , resulting in a compact representation that can be used for accurate and efficient reduced-order models and control Brunton2015amr ; Rowley2017arfm .

The majority of modal decompositions are linear taira2017aiaa ; taira2019aiaa , although emerging techniques in machine learning are providing improved nonlinear pattern extraction Brunton2020arfm . Linear regression and least-squares optimization are particularly widely used, as in the ubiquitous proper orthogonal decomposition (POD) (also known as principal component analysis (PCA)), Lumley:1970 ; aubry1988dynamics ; berkooz1993proper ; HLBR_turb ; Towne2018jfm and the emerging dynamic mode decomposition (DMD) Schmid2010jfm ; Rowley2009jfm ; Tu2014jcd ; Kutz2016book . POD provides a principled approach to decomposing high-dimensional fluid flow data into a hierarchy of orthogonal modes that are ordered in terms of their ability to capture the energy in the flow; because these modes are orthogonal, it is possible to obtain reduced-order models by Galerkin projection of the Navier-Stokes equations onto a truncated POD basis Noack2003jfm ; Carlberg2017jcp ; Loiseau2017jfm . DMD is a related technique to decompose a flow into spatiotemporal coherent structures that are each constrained to have coherent and linear dynamics in time. The least-squares regression underlying both of these approaches is highly susceptible to outliers and corrupted data candes2011jacm ; Brunton2019book , which is common in experimental measurement techniques such as particle image velocimetry (PIV). Thus, modal decomposition techniques such as POD/PCA and DMD are fragile with respect to outliers. Further, even though POD is robust to Gaussian white noise Gavish2014arxiv , DMD is sensitive to noisy data (Bagheri2014pof, ; Dawson2016ef, ; Hemati2017tcfd, ). Even techniques that are robust to corrupt velocity fields, such as finite-time Lyapunov exponents (FTLE) and Lagrangian coherent structures (LCS) (Haller2002pof, ; shadden2005pnp, ; shadden2007jfm, ; Green2011jfm, ; raben2014eif, ; Farazmand2012chaos, ), are unable to process velocity fields with missing regions, which is common in experimentally acquired data, so that interpolation must be used to fill in missing vectors. In this work, we explore the use of the robust principal component analysis (RPCA) (candes2011jacm, ) to process corrupt flow fields, leveraging global correlations in the data. We emphasize the impact of this approach on modal analysis, including POD/PCA and DMD.

1.1 Experimental challenges

Experimental techniques to measure fluid flows have evolved rapidly over the past century, with the ultimate goal of acquiring full flow fields with high spatial and temporal resolution. Laser-based imaging techniques have evolved from point measurements (yeh1964apl, ; foreman1965apl, ) to 2D and 3D field measurements (willert1991digital, ; adrian1991ar, ; adrian2005eif, ; adrian2011cup, ; westerweel2013ar, ; raffel2018particle, ). Particle Image Velocimetry (PIV) has since become a cornerstone of experimental fluid mechanics, providing non-intrusive velocity field measurements across a range of applications. Improvements in PIV hardware, including more powerful lasers, higher resolution and frame rate cameras, advanced image processing technology, and the development of tomographic PIV are providing unprecedented views into real flows. Despite the undeniable success of PIV, there are several well-known challenges to acquiring clean and accurate data. Multiple factors in the PIV data acquisition and processing pipeline can contribute to velocity vector outliers that degrade the resulting velocity fields. These include inadequate illumination and irregularities in the light field, background speckle, seeding density and non-passivity of the particles, sharp gradients in flow properties, optical issues, such as alignment and aberration, limited resolution and shot noise in the image recording, and out of plane motion of the particles when measuring in 2D (huang1997mst, ; adrian2011cup, ). Because of a fundamental tradeoff between the quantity and quality of PIV data, in both space and time, researchers continue to push the resolution limits of current systems. Thus, flow fields acquired with PIV often have missing and/or corrupt measurements. This has motivated processing techniques to improve PIV data quality of PIV data (hart2000eif, ; nogueira1997data, ; garcia2011eif, ), including predictor-corrector schemes (schrijer2008effect, ) and spatial filtering to remove frequencies not possible for the measurement resolution (discetti2013spatial, ). The identification of spurious vectors has been studied extensively (westerweel1994eif, ; duncan2010universal, ), and the normalized median filter is a robust and well-used method (westerweel2005eif, ). However, missing vectors often cluster in regions of high shear, presenting a challenge for standard vector validation and interpolation methods that rely on local flow information (westerweel2005eif, ). In this work, we leverage robust statistics and global spatiotemporal coherent structures across the entire dataset to fill in missing measurements and improve modal decomposition of fluid flow fields.

1.2 Contributions of this work

We investigate the use of RPCA (candes2011jacm, ), a robust variant of POD/PCA (Pearson:1901, ; Brunton2019book, ). RPCA uses a sparsity-promoting optimization to decompose a data matrix into the sum of a low-rank matrix containing coherent structures and a sparse matrix containing outliers. RPCA was originally popularized in the Netflix matrix completion algorithm for its recommender system (wong2015rad2, ) and has since been widely used for image and video processing (bouwmans2014cviu, ), electrical capacitance tomography (lei2013sensors, ), and voice separation (huang2012ieee, ). Here, we use RPCA filtering to process flow field data from several simulations and experiments. Figure 1 demonstrates the ability of RPCA to uncover and isolate the dominant low-rank coherent structures from sparse outliers in flow data from an idealized example. In addition to directly analyzing and processing flow field data, we also perform PCA and DMD modal analyses on the data before and after RPCA filtering to assess its performance.

Here we consider a range of simulated and experimentally acquired flow fields of varying complexity to isolate and analyze various aspects of the algorithm applied to data from fluid mechanics. First, we investigate high-fidelity flow fields from direct numerical simulations of a laminar flow past a cylinder and a turbulent channel flow, where it is possible to artificially add corrupt velocity field vectors to compare the RPCA filtered fields with a known ground truth. Next, we apply the method to two experimentally acquired datasets, including a companion laminar flow past a cylinder and measurements of a cross-flow turbine wake. Although there is not a ground truth model for these flows, it is possible to estimate the effect of RPCA filtering in reducing outliers and corruption by analyzing the DMD spectrum, which has well-stereotyped behavior for such periodic wake flows Bagheri2014pof . In all cases, we show that the RPCA filtered fields yield DMD spectra that are more consistent with a periodic wake in the absence of noise.

This work is organized as follows: First, we present the standard POD/PCA and DMD modal analysis techniques in Sec. 2, followed by the RPCA method in Sec. 3. Section 4 describes the four flow fields used in this analysis. Results of RPCA filtering on these flow fields and its impact on downstream modal analysis are presented in Sec. 5.

2 Modal analysis

Extracting coherent structures from high-dimensional data has been a central challenge in fluid mechanics for decades. Here we review two leading modal decomposition techniques for data from fluid mechanics, the proper orthogonal decomposition (POD), also known as principal component analysis (PCA) (Sec. 2.1), and the dynamic mode decomposition (DMD) (Sec. 2.2). Both methods apply equally well to data from simulations or experiments. We use these two modal decomposition techniques to assess the effectiveness of RPCA filtering, in processing and correcting corrupt flow fields. Both techniques are based on the singular value decomposition (SVD) Golub1980siamjna ; Golub1970nm ; Businger1969cacm ; Golub1965siamb and there are several detailed discussions of these decompositions Kutz2016book ; taira2017aiaa ; taira2019aiaa ; Brunton2019book . The RPCA algorithm is explained in Sec. 3.

In this work, we follow the flowchart shown in Fig.  2. RPCA filtering is applied to a data matrix (X) which is bisected into the low rank structure (L) and sparse (S) subspaces. From there, RPCA modes and DMD modes are calculated from the low rank data. We also calculate the POD/PCA and DMD modes on the data matrix.

2.1 Proper orthogonal decomposition (POD)

Proper orthogonal decomposition–referred to as principal component analysis (PCA) throughout results—is a widely used method to identify spatially correlated coherent structures from data, decomposing the flow field into a linear combination of orthogonal modes that are arranged hierarchically by energy content. There are several variants of POD taira2017aiaa ; Lumley:1970 ; Sirovich:1987 ; berkooz1993proper ; HLBR_turb ; Towne2018jfm , and we will present a variant of the snapshot POD of Sirovich that relies on the numerically stable SVD Sirovich:1987 ; Brunton2019book . First, flow field data (e.g., a velocity or vorticity field) is measured or computed on a discrete spatial grid, and $m$ snapshots of these flow fields are collected at various times $t_{1},t_{2},\cdots,t_{m}$. The flow field data at time $t_{k}$ may be reshaped into a column vector ${\bf x}_{k}={\bf x}(t_{k})\in\mathbb{R}^{n}$, where $n$ denotes the number of flow variables times the number of spatial grid locations. Next, a data matrix is formed by arranging the column vectors ${\bf x}_{k}$ in a matrix ${\bf X}$:

[TABLE]

Finally, POD modes are obtained by computing the singular value decomposition of ${{\bf X}\in\mathbb{R}^{n\times m}}$:

[TABLE]

where T defines the matrix transpose, ${\bf U}\in{\mathbb{R}}^{n\times n}$, ${\bf\Sigma}\in{\mathbb{R}}^{n\times m}$ and ${\bf V}\in{\mathbb{R}}^{m\times m}$. The columns of ${\bf U}$ are POD modes with the same dimension as a flow field ${\bf x}$. POD modes are orthonormal so that ${\bf U}^{T}{\bf U}={\bf I}$; similarly ${\bf V}^{T}{\bf V}={\bf I}$. Moreover, the columns of ${\bf U}$ (resp. rows of ${\bf V}^{T}$) are arranged in order of their importance in describing the data. The importance of each mode (i.e., column of ${\bf U}$) is given by the corresponding entry of the non-negative, diagonal matrix of singular values $\boldsymbol{\Sigma}\in\mathbb{R}^{n\times m}$.

The matrix ${\bf X}$ will exhibit low-rank structure, so that it is well approximated by the first $r\ll m<n$ columns of ${\bf U}$ and ${\bf V}$:

[TABLE]

where ${\bf U}_{r}$ and ${\bf V}_{r}$ denote the first $r$ columns of each matrix and $\boldsymbol{\Sigma}_{r}$ denotes the first $r\times r$ sub-block of $\boldsymbol{\Sigma}$. In fact, the Eckart-Young theorem states that this is the optimal rank-$r$ approximation of the matrix ${\bf X}$ in a least-squares sense. More details about the SVD can be found in Brunton2019book .

After truncating all but the first $r$ dominant modes, a flow field snapshot ${\bf x}$ may be approximated by a linear combination of these modes:

[TABLE]

where $\alpha_{k}$ is the POD mode coefficient.

Because the POD modal basis is orthogonal, it is possible to obtain a reduced-order nonlinear dynamical system for the evolution of the coefficients $\alpha_{k}(t)$ in time via Galerkin projection of the Navier-Stokes equations onto the POD basis. In this way, the POD basis may be thought of a data-driven generalization of the Fourier basis that is tailored to a particular flow field. POD is also closely related to principal component analysis (PCA) Pearson:1901 , the Karhunen–Loève decomposition karhunen1947lineare , empirical orthogonal functions lorenzMITTR56 , or the Hotelling transform hotellingJEdPsy33_1 ; hotellingJEdPsy33_2 .

2.2 Dynamic mode decomposition

DMD is a modal decomposition technique that simultaneously identifies spatially coherent modes that are constrained to have the same linear behavior in time, given by oscillations at a fixed frequency with growth or decay Schmid2010jfm ; Kutz2016book . Thus, the dynamic mode decomposition provides a dimensionality reduction into a set of spatial modes along with a linear model for how these modes evolve in time. This is in contrast to POD, which results in orthogonal modes arranged in terms of energy content and without consideration of dynamics. However, in many formulations, DMD is closely related to POD, and may be thought of as a linear combination of POD modes that results in linear evolution in time. DMD also has deep connections to nonlinear dynamical systems via Koopman operator theory Mezic2005nd ; Rowley2009jfm ; Tu2014jcd ; Mezic2013arfm ; Kutz2016book .

In the original formulation of DMD Schmid2010jfm , the snapshots in the data matrix in Eq. (1) are spaced evenly in time, so that $t_{k}=k\Delta t$ with $\Delta t$ sufficiently small to resolve the highest frequencies in the dynamics. Generalizations have since been formulated to allow for non-sequential time-series Tu2014jcd ; Askham2018siads and for data that is under-resolved in space Brunton2015jcd ; Gueniat2015pof or time Tu2014eif ; however, for simplicity, we will present the standard exact DMD formulation of Tu et al. Tu2014jcd with evenly spaced and sequential snapshots. DMD seeks to identify the leading eigenvalues and eigenvectors of the best-fit linear operator ${\bf A}$ that evolves snapshots forward in time:

[TABLE]

The eigenvectors $\boldsymbol{\phi}$ of $\mathbf{A}$ have the dimensions of a flow-field and correspond to spatiotemporal coherent structures whose dynamics in time evolve according to the associated eigenvalue $\gamma$.

In practice, this operator is identified by first splitting the data in Eq. (1) into two matrices:

[TABLE]

and then solving for the best fit operator that satisfies

[TABLE]

via the following least-squares optimization problem

[TABLE]

Here we are minimizing the Frobenius norm $\|\cdot\|_{F}$ via the pseudo-inverse ${{\bf X}^{\dagger}={\bf V}\boldsymbol{\Sigma}^{-1}{\bf U}^{T}\approx{\bf V}_{r}\boldsymbol{\Sigma}_{r}^{-1}{\bf U}_{r}^{T}}$.

In practice, the matrix ${\bf A}$ is far too large to analyze directly, and instead, we project ${\bf A}$ onto an $r$-dimensional POD subspace, given by the columns of ${\bf U}_{r}$:

[TABLE]

The eigenvalues of $\tilde{\bf A}$ are the same as the eigenvalues of ${\bf A}$, which are known as the DMD eigenvalues. They are computed via the eigendecomposition of the $r\times r$ matrix $\tilde{\bf A}$:

[TABLE]

Finally, the corresponding DMD modes are reconstructed using the full-dimensional data along with the reduced eigenvectors in ${\bf W}$:

[TABLE]

This formula for the eigenvectors is from the exact DMD algorithm Tu2014jcd ; Kutz2016book ; the original formulation of Schmid Schmid2010jfm computes modes as $\boldsymbol{\Phi}={\bf U}_{r}{\bf W}$.

With the DMD modes $\boldsymbol{\Phi}$ and eigenvalues $\boldsymbol{\Gamma}$ it is possible to reconstruct the state at time $k\Delta t$

[TABLE]

where the vector ${\bf b}$ of mode amplitudes is generally computed as

[TABLE]

More principled approaches to select the few dominant modes have been considered based on sparsity-promoting optimization Jovanovic2014pof .

The spectral expansion above may also be written in continuous-time by introducing the continuous eigenvalues $\omega=\log(\gamma)/\Delta t$:

[TABLE]

where $\boldsymbol{\Omega}$ is a diagonal matrix containing the continuous-time eigenvalues $\omega_{j}$.

DMD is known to be extremely sensitive to noisy data (Bagheri2014pof, ; Dawson2016ef, ; Hemati2017tcfd, ), and the eigenvalues specifically suffer from a bias that is not reduced with increasing data. There are several modifications to make DMD more robust to noise, including averaging forward-time and backward-time operators Dawson2016ef , total least squares Hemati2017tcfd , and variable projection Askham2018siads . For periodic wake data, as explored in three of the examples in this paper, the discrete-time eigenvalues should occur in complex conjugate pairs $\gamma,\bar{\gamma}$ exactly on the unit circle in the complex plane for clean data. Similarly, the continuous-time eigenvalues should be in complex conjugate pairs $\pm i\omega$ on the imaginary axis, indicating pure oscillations with no growth or decay Bagheri2013jfm ; Bagheri2014pof . In Bagheri2014pof , Bagheri characterized the perturbative effect of noise on these eigenvalues, deriving an asymptotic expression for how high frequency eigenvalues become increasingly affected by noise. If the true continuous-time eigenvalue should be $\pm i\omega$ in the absence of noise, Bagheri showed that in the presence of perturbative white noise with magnitude $\epsilon\ll 1$ the observed eigenvalue pair is

[TABLE]

where $C$ is a sensitivity constant. Thus, low noise levels cause a spurious real-valued damping $-\epsilon C\omega^{2}$ that is quadratic in the frequency. We will make extensive use of this property to assess the quality of our RPCA filtered fields by computing the ratio of the best-fit factor $\epsilon C$ before and after applying RPCA filtering. $\epsilon C$ should decrease as a consequence of reduced noise and corruption.

3 Robust extraction of fluid coherent structures

Techniques based on least-squares regression, such as POD/PCA and DMD, are highly susceptible to outliers and corrupted data, making them fragile with respect to some experimental measurement errors. To mitigate this sensitivity, Candes et al. candes2011jacm have developed a robust principal component analysis (RPCA) that seeks to decompose a data matrix ${\bf X}$ into a structured low-rank matrix ${\bf L}$ that is characterized by dominant coherent structures and a sparse matrix ${\bf S}$ containing outliers and corrupt data:

[TABLE]

The principal components of ${\bf L}$ are robust to outliers and corrupt data, which are isolated in ${\bf S}$. This decomposition, also referred to as a filter, has profound implications for many modern problems of interest, including video surveillance (where the background objects appear in ${\bf L}$ and foreground objects appear in ${\bf S}$), facial recognition (eigenfaces are in ${\bf L}$ and shadows, occlusions, etc. are in ${\bf S}$), natural language processing and latent semantic indexing, and ranking problems111The ranking problem may be thought of in terms of the Netflix prize for matrix completion. In the Netflix prize, a large matrix of preferences is constructed, with rows corresponding to users and columns corresponding to movies. This matrix is sparse, as most users only rate a handful of movies. The Netflix prize seeks to accurately fill in the missing entries of the matrix, revealing the likely user rating for movies the user has not seen..

Mathematically, the goal is to find ${\bf L}$ and ${\bf S}$ that satisfy the following:

[TABLE]

However, neither the $\text{rank}({\bf L})$ nor the $\|{\bf S}\|_{0}$ terms are convex, making this optimization intractable. Similar to the compressed sensing problem, it is possible to solve for ${\bf L}$ and ${\bf S}$ with high probability using a convex relaxation of (16):

[TABLE]

where $\|\cdot\|_{*}$ is the nuclear norm, given by the sum of singular values which is a proxy for the rank of the matrix and ${\lambda_{0}=\lambda/\sqrt{\max(n,m)}}$ and $\|\cdot\|_{1}$ is the 1-norm of the matrix. The solution to (17) converges to the solution of (16) with high probability if $\lambda=1$, where $n$ and $m$ are the dimensions of ${\bf X}$, given that ${\bf L}$ is not sparse and ${\bf S}$ is not low-rank. In the examples below, these assumptions may only be partially valid, so the optimal value of $\lambda$ may vary slightly. The convex problem in (17) is known as principal component pursuit (PCP), and may be solved using the augmented Lagrange multiplier (ALM) algorithm. Specifically, an augmented Lagrangian may be constructed:

[TABLE]

Where $Y$ is the matrix of Lagrange multipliers and $\upsilon$ is a hyperparameter. We then solve for ${\bf L}_{k}$ and ${\bf S}_{k}$ to minimize $\mathcal{L}$, update the Lagrange multipliers

[TABLE]

and iterate until convergence. In this work, an inexact ALM implementation from sobral2015lrs is used. The alternating directions method (ADM) (Lin2010adm, ; Yuan2009adm, ) provides another simple procedure.

After the low-rank matrix ${\bf L}$ is obtained, it is possible to compute robust POD/PCA modes (Sec. 2.1) as in Eq. (2):

[TABLE]

Henceforth, we refer to the modes in ${\bf U}$ from ${\bf UL}$ as RPCA modes. Similarly, it is also possible to compute robust DMD (RDMD) modes and eigenvalues.

4 Model flows

We demonstrate RPCA filtering on several example data sets of varying complexity, drawn from direct numerical simulations (DNS) and PIV data from experiments. Figure 3 provides an overview of the four example flow fields.

4.1 Cylinder flow

Flow past a cylinder is a canonical example in fluid mechanics. We consider data from DNS at a diameter-based Reynolds number of 100 and from PIV measurements at Reynolds number 413 (Shang2015flexibility, ).

The DNS data is generated by simulating the two-dimensional incompressible Navier-Stokes equations using the immersed boundary projection method taira:fastIBPM ; taira:07ibfs . The computational domain comprises four nested grids: the finest grid covers a domain of $9\times 4$ and the largest grid covers a domain of $72\times 32$, where lengths are non-dimensionalized by the cylinder diameter. Each grid contains $449\times 199$ points with a resolution of 50 points per cylinder diameter. The time-step is $\Delta t=0.02$ and data is sampled at intervals of $10\Delta t$ (30 times the vortex shedding frequency) with $m=150$ snapshots saved, covering $5$ vortex shedding cycles. The DNS provides a benchmark, where the uncorrupted flow field is known, to quantitatively assess performance of RPCA filtering on data with artificial salt & pepper corruption. Corrupted sample points are chosen uniformly in space and time at a given rate, and both the $u$ and $v$ velocity components at each selected location are randomly assigned a value of $\pm 10$ times the standard deviation of the streamwise velocity data. In addition, we consider a second case where corrupted sample points are chosen with a bias towards regions of high vorticity or shear magnitude, which is more physically realistic for PIV data. In this second case, we select measurements for corruption based on a probability density given by $\alpha+|\boldsymbol{\omega}|$, where $\alpha$ is a small positive constant and $|\boldsymbol{\omega}|$ is the absolute value of the vorticity; when the rate of corruption is sufficiently high, these corrupted fields begin to resemble the uniformly corrupted cases, but with more corruption in vortex cores. Because vorticity is calculated from velocity fields using a finite-difference derivative, there is a higher rate of corruption in the vorticity fields than in the velocity fields.

The PIV data has frame size of $135\times 80$ grid points with a resolution of $8$ points per cylinder diameter. Data is sampled at a rate of $20\,$Hz (125 times the shedding frequency) with $m=8,000$ snapshots saved, which corresponds to $64$ vortex shedding cycles.

4.2 Turbulent channel flow

For a more complex and multi-scale flow, we consider DNS data from a forced, fully developed turbulent channel flow data with a friction velocity Reynolds number of $Re_{\tau}=1,000$, from the Johns Hopkins Turbulence Database (JHTDS) (graham2016jot, ). This example provides a test case to see how turbulent kinetic energy at various scales is filtered depending on the level of added noise. The addition of noise is similar to the cylinder DNS where randomly selected sample points of the streamwise and cross-stream velocity fields are assigned a value of $\pm 10$ standard deviations of the streamwise velocity data. Due to the size of the full dataset, we only consider two-dimensional fields on the mid-plane, with a $512\times 512$ grid of three component velocity measurements spanning the channel width. Data is sampled at a rate of 966 times the mean flow-through time with $m=1000$ snapshots.

4.3 Cross-flow turbine wake

Finally, we consider PIV wake data from a cross-flow turbine experiment conducted at the University of Washington. Cross-flow turbines can be used to extract power from wind and water currents for renewable energy generation. This flow exhibits both coherent and broadband phenomena, and provides a challenging test-case RPCA filtering. The frame consists of $158\times 98$ grid points, at a resolution of 99 points per rotor diameter. Data is sampled at a rate of 32 times the blade-pass frequency with $m=1000$ snapshots. Vectors were calculated using a multi-grid, multi-pass algorithm with adaptive image deformation (scarano2001mst, ). Resulting vector fields were then validated using a normalized median filter with potential replacement by secondary correlation peaks. The cross-correlation and validation steps result in missing data, particularly in regions of high vorticity and shear. To apply RPCA filtering, these missing values are randomly assigned a value of $\pm 10$ standard deviations of the streamwise flow data, in contrast to the experimental cylinder wake where missing measurements were previously interpolated.

5 Results

We now explore the ability of RPCA filtering to isolate and remove noise and corruption from the example flow fields. We will begin with the simulated and experimental flow past a cylinder, followed by data from the Johns Hopkins turbulent channel flow simulation, and ending with the experimental wake of a cross-flow turbine.

5.1 Cylinder flow

Figure 4 shows the results of RPCA filtering for flow past a cylinder, providing a side-by-side comparison of PIV and corrupted DNS data data. Although the Reynolds numbers differ by a factor of four, the flow fields are qualitatively similar, characterized by periodic, laminar vortex shedding. For $\lambda=1$, the data is correctly segmented with the coherent flow in ${\bf L}$ and the sparse corruption in ${\bf S}$. When $\lambda$ is too small, RPCA filtering is overly aggressive, incorrectly including relevant flow structures in ${\bf S}$, and when $\lambda$ is too large, the corruption is not filtered.

For the experimental data in Fig. 4 (right), the optimal value of $\lambda$ is less clear. For $\lambda=0.1$, the low-rank field ${\bf L}$ is visually smoother than the field at $\lambda=1$, but the sparse matrix ${\bf S}$ contains a significant portion of the wake structures, indicating over-filtering. This filtering becomes more pronounced in the movies, where it is clear that much of the high-frequency “noise” in the bypass flow is actually free-stream turbulence, which is consistent with the turbulence intensity of the experiments. Further, as subsequently discussed, when we compute the RPCA modes, it is clear that the $\lambda=0.1$ case is heavily filtering out all but the first three modes. Thus, it appears that the theoretically optimal value $\lambda=1$ has the best performance, although there may be a tradeoff between filtering ambient free-stream turbulence and coherent structures of interest in experiments.

Figure 5 shows the results of RPCA filtering on the simulated data for the second case of vorticity-biased corruption. Again, in all cases, the theoretically optimal value of $\lambda=1$ yields the best segmentation of the corruption into the matrix ${\bf S}$. When the rate of corruption is increased from $1\%$ (left) to $10\%$ (right), the free-stream flow begins to become corrupted, resembling the uniform corruption case in Fig. 4. The mean error and relative nuclear norm of the low-rank matrix (${\bf L}$) compared to the true, uncorrupted data (${\bf X}$) are shown in Fig. 6 for varying percentages of corrupt entries. Statistically, results are similar for vorticity-biased and randomly-distributed corruption. In both cases, RPCA filtering is remarkably robust to corruption, even for corruption in excess of $50\%$ of the measurements. This laminar vortex shedding example is an ideal application for RPCA filtering, as the true flow field is low rank and the corruption is sparse; it is unlikely that this will hold as well for data exhibiting broadband turbulence.

5.1.1 PCA analysis for cylinder wake flows

We now investigate the impact of RPCA filtering on modal decompositions. Figures 7 and 8 show four leading PCA and RPCA modes for $1\%$ and $10\%$ vorticity-biased corruption, respectively. The first mode corresponds to the mean flow, and the remaining modes come in energetic pairs where the corresponding coefficients $\alpha_{2j}$ and $\alpha_{2j+1}$ oscillate sinusoidally at the same frequency but $\pi/2$ out of phase, sweeping out a circle in the phase plane. Thus, we only show one mode, $\boldsymbol{u}_{2j+1}$ from each of the first three energetic mode pairs. In all cases, the RPCA modes show dramatic improvement, while significant artifacts remain in the PCA modes. We also investigate the effect of increasing the amount of data, and there is a clear improvement in RPCA modes from 2 to 5 vortex shedding cycles; in contrast, the PCA modes do not improve appreciably with more data.

To quantify the improvement observed above, we compute the $L_{2}$ error between the PCA and RPCA modes of corrupted data and the PCA modes for the clean data (i.e., DNS results) as a function of the number of shedding periods. As show in Fig. 9, the RPCA mode velocity fields quickly converge to a small error as the amount of data is increased for both the $1\%$ and $10\%$ corruption cases, while the PCA modes converge much more slowly and still have considerable error after 5 shedding periods are included in the analysis.

The modes for the PIV data for the cylinder flow are shown in Fig. 10. This figure highlights the effect of $\lambda$, the sparsity hyperparameter, which was previously discussed with respect to Fig. 4. In this case, we do not have a clean ground-truth data set to compare against. Although the flow field in Fig. 4 appears to have less corruption for $\lambda=0.1$, here we see that all RPCA modes after the first three modes are heavily filtered, as seen in the rapid drop off in the singular values after the third mode. The corresponding modes are highly corrupt, further supporting that $\lambda=0.1$ is not a good choice. In contrast, the RPCA modes for the theoretically optimal $\lambda=1$ case appear to have slightly less free-stream corruption than the PCA modes. Also, as expected, for a large enough value of $\lambda$, the RPCA filtering has little effect on the modes.

5.1.2 DMD analysis for cylinder wake flows

DMD is known to be quite sensitive to noisy data, making this a challenging test case for RPCA filtering. Figure 11 shows the discrete-time eigenvalues for the cylinder DNS data with vorticity-biased corruption. For the cylinder wake, the uncorrupted or true DMD eigenvalues may be computed from the noiseless data, and they are equally spaced on the unit circle in the complex plane. In all cases, the RPCA-filtered DMD (RDMD) data results in dramatically better agreement with the uncorrupted or true DMD eigenvalues compared with the corrupted DMD eigenvalues. Even with only a single period of data and $\eta=10\%$ corruption, the RDMD values capture the first six low-frequency mode pairs; in contrast, even with five periods of data and as little as $\eta=1\%$ corruption, corrupted DMD only captures the first two low-frequency mode pairs, and with considerably more spurious damping. To see this more clearly, we plot the eigenvalues in continuous-time in Fig. 12, where the $x$-axis is the imaginary eigenvalue component and the $y$-axis is the real eigenvalue component, which is a standard way to plot DMD eigenvalues Schmid2010jfm . Here, the best-fit parabolas for the RDMD eigenvalues and the first seven corrupted DMD eigenvalues are shown in dashed lines. The curvature of these parabolas is directly related to the noise amplitude, as in Eq. (14) from Bagheri2014pof . The same continuous-time eigenvalue plot is shown for the PIV cylinder wake data in Fig. 13. In both cases, we see that the parabolic eigenvalue fit for the RDMD eigenvalues has a smaller curvature than for the corrupted DMD eigenvalues, indicating a quantitative and significant reduction in noise.

5.2 Turbulent channel flow

The tradeoff between filtering corruption and small-scale structures is also apparent in the turbulent channel flow DNS. Unlike the cylinder wake, this flow field contains broadband turbulent phenomena across multiple spatial and temporal scales. Figure 14 shows RPCA filtering for various levels of corruption, sweeping across the tuning parameter $\lambda$. The corresponding turbulent kinetic energy (TKE) is shown in Fig. 15, providing a summary of the various scales that are filtered. The value of $\lambda$ that preserves the true TKE spectrum varies with the degree of velocity field corruption. In the uncorrupted case ($\eta=0$), we can clearly see the effect of filtering on the turbulent coherent structures, indicating that some fine-scale structures are filtered for $\lambda=1$. As the degree of corruption increases to $\eta=2\%$, we see that the curves for $\lambda\leq 2$ remain relatively unchanged, although the $\lambda>2$ curves begin to exhibit spurious high-frequency spatial structures (i.e., corruption is present in $\bf L$). As the rate of corruption increases to $\eta=10\%$, spurious high-frequency energy also appears for $\lambda=2$. In this case, it is clear that the optimal filtering value $\lambda$ changes with the level of corruption. For relatively limited corruption, a larger value of $\lambda$ may be used, but must be decreased towards the theoretically optimal value of $\lambda=1$ for higher levels of corruption. Finally, we note that, unlike the cylinder wake cases, it is not surprising that $\lambda=1$ is sub-optimal because the channel flow is not fundamentally low rank which deviates from an underlying assumption of the RPCA algorithm.

5.3 Cross-flow turbine wake

As a final example, we consider the use of RPCA filtering to identify outliers and fill in missing PIV data collected in the wake of a cross-flow turbine, as shown in Fig. 16. There are several stages in the PIV processing pipeline where RPCA filtering could be applied, including after initial cross-correlation, after conventional normalized median filter vector validation, and after linear interpolation. For the cases shown here, we use $\lambda=1.6$, which results in a velocity in the bypass flow, or lower third of the frame, that visually matches the frequency content of the unfiltered data. There are enough missing velocity vectors (23% and 20%, respectively) to degrade the effectiveness of both median filtering and interpolation. In contrast, RPCA filtering produces flow fields that capture dominant coherent structures for either cross-correlated or median-filtered fields. Finally, by investigating the standard deviation of all flow fields collected at a given turbine angular position (i.e., phase), it is clear that the RPCA filtering can be used to remove artifacts introduced by linear interpolation. This is consistent with the intuition that vector validation and interpolation should fail in these regions where there is high density of missing data is spatially clustered.

The continuous-time DMD eigenvalues for the cross-flow turbine wake are shown in Fig. 17. In this plot, the parabolic fits for DMD eigenvalues computed after interpolation and RDMD eigenvalues computed after vector validation are displayed as dotted lines. For this case, the coefficient of the parabolic fit for the RDMD-based eigenvalues is six times smaller than the parabolic fit coefficient for the DMD-based eigenvalues. This demonstrates a significant quantitative improvement of the DMD spectrum using RPCA filtering to process the data.

6 Conclusions and discussion

In this work, we have demonstrated the ability of RPCA filtering to effectively recover dominant coherent structures from corrupt flow fields with missing measurements. Unlike standard POD/PCA, which is based on least squares and is susceptible to outliers and corruption, RPCA utilizes sparse optimization to separate a data matrix into a low-rank matrix containing correlated structures and a sparse matrix containing the spurious entries.

We apply RPCA filtering to several types of fluid flow data (DNS and PIV), ranging from laminar vortex shedding behind a circular cylinder, to fully turbulent channel flow DNS, and concluding with an experimental flow past a cross-flow turbine. These flows exhibit a variety of phenomena and a range of measurement quality. The DNS examples provide us with a baseline, where it is possible to add corruption to quantitatively assess the performance of RPCA. For flow past a cylinder in DNS, RPCA filtering is extremely effective at separating the true flow field from considerable corruption, with robust recovery even in flow fields with excess of 50% of the measurements corrupted. In the experimental counterpart, RPCA is still able to remove large outliers and corruption, although there is a tradeoff between filtering the background turbulence and coherent structures in the wake. The fully turbulent channel flow DNS provides an opportunity to more fully explore this tradeoff in a controlled setting, where we can incrementally increase the corruption ratio and observe the filtering effects on various spatial frequencies. As expected, an increasingly aggressive filtering leads to degradation at higher wavenumbers, although dominant coherent structures are robustly preserved. Finally, the wake behind a cross-flow turbine provides a practical real-world flow that directly benefits from improved PIV processing. In all three wake flows we also assess the performance of RPCA filtering to yield more accurate modal decompositions. Although we do not have ground truth measurements and modal decompositions, except in the case of direct numerical simulations, we know that continuous-time DMD eigenvalues should be arranged on the imaginary axis in the complex plane for clean data, and deviations from this may be quantified using the derivation from Bagheri Bagheri2014pof . In all three cases, we see considerable reduction in spurious damping, indicating the de-noising effectivness of RPCA. Based on these results, we believe that RPCA can be a valuable algorithm in the arsenal of PIV processing and filtering techniques, particularly when the processing pipeline culminates in modal analysis.

There are a number of future directions motivated by this work. First, RPCA depends on the hyperparameter $\lambda$, and a better understanding of how to objectively choose $\lambda$ for different conditions is important. Because RPCA is based on sparse, non-convex optimization, it is also likely that improved optimization techniques may improve speed and robustness. Although this work considered three-dimensional flows, the data comprised two-dimensional cross-sections, and the current analysis could be extended to flow volumes. In principle, the RPCA method should generalize, although there may be computational scaling challenges. It would also be useful to extend this work to PIV measurements of other turbine configurations (posa2016wake, ). Finally, the quality of the RPCA filtered flow fields for additional downstream analyses should be assessed for example, in dynamical systems modeling via Galerkin projection (Noack2003jfm, ) or regression (Loiseau2017jfm, ) onto the filtered modes and in control (berger2015active, ).

Acknowledgments

We gratefully acknowledge funding from the Army Research Office (ARO W911NF-19-1-0045) and Naval Facilities Engineering Command. We would also like to thank Jared Callaham, Nathan Kutz, Kazuki Maeda, and Joshua Proctor for valuable discussions.