Physics-aware nonparametric regression models for Earth data analysis

We used the SeaBAM dataset [30, 31], which gathers 919 in-situ measurements of chlorophyll concentration around the United States and Europe. The dataset contains in situ pigments and remote sensing reflectance measurements (Rrs, [sr⁻¹]) at a set of given wavelengths (412, 443, 490, 510 and 555 nm) that are present in the SeaWiFS ocean color satellite sensor. The chlorophyll concentration values are between 0.019 and 32.79 mg m⁻³. Although SeaBAM data originate from various researchers, the variability in the radiometric data is limited. At high Chla concentrations CC (mg m⁻³), the dispersion of radiance ratios increases, mostly because of the presence of more optically complex (Case II) waters, that is low values of the ratio of pigment concentration to scattering coefficient. At lowest concentrations the highest $R(490)/R(555)$ ratios are slightly lower than the theoretical limit for clear natural waters. More information about the data can be obtained from SEABAM, and an extensive analysis in [31]. In addition to the observational data, we used several models to guide PKL. Morel1 and CalCOFI 2-band linear are described by $s = 10^{a_0+a_1\kappa}$, while OC2/OC4 models follow $s = a_0 + 10^{a_1+a_2\kappa+a_3\kappa^2+a_4\kappa^3}$. The ratio κ depends on the physical model used [31]: the Morel1 model uses $\kappa = \log(R(443)/R(555))$, the CalCOFI and OC2 models use $\kappa = \log(R(490)/R(555))$, while the OC4 model uses $\kappa = \log(\max\{R(443),$ $R(490),R(510)\}/R(555))$. The goal here is to predict the concentrations from reflectance measurements while being consistent with these parametric models that encapsulate some domain knowledge.

We collected hyperspectral data from the CHRIS sensor as well as in situ measurements of chlorophyll content (Chla), leaf area index (LAI) and fractional vegetation cover (fCOVER) [32]. LAI is a dimensionless quantity that characterizes plant canopies. It is defined as the one-sided green leaf area per unit ground surface area in broadleaf canopies. fCOVER corresponds to the fraction of ground covered by green vegetation, and in practice it quantifies the spatial extent of the vegetation. It is indeed independent from the illumination direction and sensitive to the vegetation amount and derived from the LAI plus a number of structural parameters of the canopy. The data were obtained during the SPARC-2003 (SPectra bARrax Campaign) and SPARC-2004 campaigns in Barrax, Spain. The region consists of approximately 65% dry land and 35% irrigated land. Green LAI was derived from canopy measurements made with a LiCor LAI-2000 digital analyzer. Each elementary sampling unit (ESU) was assigned to a LAI value, which was obtained by the average of 24 measures (8 data readings × 3 replications) [33]. fCOVER was estimated from ground measurements using hemispherical photographs taken with a digital camera with a fish-eye lens. The final fCOVER estimate for each ESU was calculated as the average of twelve measurements. In total, nine crop types (garlic, alfalfa, onion, sunflower, corn, potato, sugar beet, vineyard and wheat) were sampled, with field-measured values of LAI that vary between 0.4 and 6.3, Chla between 2 and 55 µg cm⁻² and fCOVER between 0 and 1. Additionally, 30 random bare soil spectra with a biophysical (Chla, LAI, fCOVER) value of zero were added to broaden the dataset to non-vegetated samples. Concurrently, we used CHRIS images Mode 1 (62 spectral bands, 34 m spatial resolution at nadir). The images were geometrically and atmospherically corrected. A total of n = 136 data points in a 62-dimensional space were thus used to fit a PKL model. The goal here was to estimate LAI while being consistent with fCOVER, and estimate chlorophyll content while being consistent to LAI.

The US CLIVAR Working Group on large ensembles (LEs) contains a data archive of initial-condition LEs conducted with different climate models run within their CMIP5 setup. From the seven climate models available, we selected three different ones: CanESM2 for training, CSIRO-Mk3-6-0 for validation, and CESM1-CAM5 for testing. We emphasize that different train/validation/test splits are possible (and should be tested in real-world applications), but here we only show an illustration which is why we present this setup and given that results are robust for alternative choices. Model simulations for all three contained the period between 1950 and 2100 with historical and RCP 8.5 forcing conditions and incorporated more than 30 different ensemble members [34]. All model runs were defined with a spatial regridded world map of $5^\circ \times 5^\circ$ resolution (in total d = 2592 grid points). The aggregated ensemble of each model, and all the grid points, defined the predictors matrix. We split it into a training, validation and test sets. We selected the spatially explicit simulated monthly near-surface air temperature (TAS) as our variable of interest (i.e. predictors). The ENSO internal variability is captured by the Niño3.4 index, and is taken from the climate variability diagnostics package for LEs developed by NCAR's Climate Analysis section [35]. Our goal is to predict the forced climate response from a single ensemble member, and the 'true' forced response was extracted as the average across the full ensemble, which is a standard procedure in the field. The predictors and metrics used were taken as the DJF seasonal average and standardized to zero mean and unit standard deviation.

A classical problem in remote sensing and geosciences involves estimation of biophysical parameters of interest from remote sensing (satellite) observations. We are given multidimensional measurements ${\mathbf{x}}\in\mathbb{R}^d$, acquired from a satellite sensor at d different wavelengths, from which we want to estimate a parameter of interest $y\in\mathbb{R}$. This can be done using satellite and in-situ measurement pairs $({\mathbf{x}},y)$. Traditionally, the community has designed a great many empirical, parametric models for estimation. We aim to use the PKL method to solve the fitting problem and, more importantly, to assess the most appropriate parametric model to rely on by studying the consistency-vs-accuracy path diagrams.

We first illustrate the performance of the PKL model for the estimation of in-situ chlorophyll concentration from multispectral reflectance images. The two measurements are subject to high levels of uncertainty, owing both to the difficulties in ground-truth data acquisition, and the noise inherent in satellite-obtained data. Moreover, there is commonly a time mismatch between the acquired image and the recorded in-situ measurements, which is critical, for instance, for coastal water monitoring. We use the SeaBAM dataset [30, 31] (see details in Data section), and apply PKL to model ocean chlorophyll concentration y from radiances ${\mathbf{x}}: = [R(\lambda_1),\ldots,R(\lambda_d)]\in \mathbb{R}^d$ at a set of given wavelengths, λ_i , while encouraging consistency with four standard parametric models (Morel1, CalCOFI 2-band linear, OC2 and OC4) [31] in a set of four different experiments, one for each model as the ancillary variable $\mathsf{s}$. Results are shown in figure 1. PKL solves the fitting problem with improved accuracy, i.e. a lower error (RMSE) can be achieved compared to standard kernel ridge regression (µ = 0), and allows to quantify the quality of the different parametric models considered through their respective consistency-vs-accuracy paths. Encouraging consistency with more recent models, such as OC2 and OC4, leads to a lower PKL prediction error (lower RMSE) compared to older models, such as Morel1 and CalCOFI, meaning that they fit the data better (i.e. the consistency between the data and the model is higher). When increasing dependence with each parametric model, the HSIC regularizer begins to dominate the training cost function and similarity of the target variable with the ancillary variable is reinforced. This happens until a certain turning point in µ where too much dependency leads to an increased error. In summary, this simple example shows how PKL may benefit the estimation of Earth system parameters such as the in-situ chlorophyll concentration from remote sensing by enforcing consistency with well-known parametric models that encapsulate domain knowledge.

Zoom In Zoom Out Reset image size

Very often one does not have access to a parametric model as in the previous example, but to a set of ancillary observational data that the predictions should be consistent with. This is a standard case in remote sensing where some variables are coupled, e.g. greenness and chlorophyll content of the vegetation. We illustrate the use of the PKL in two cases: (1) to estimate leaf area index ($\mathsf{y} : =$ LAI) while being consistent with the fraction of green vegetation cover (${\mathsf{s}}: =$ fCOVER), and (2) to predict chlorophyll-a concentration ($\mathsf{y}: =$ Chla), while being consistent with the leaf area index (${\mathsf{s}}: =$ LAI). We use data from a terrestrial campaign where hyperspectral images and the aforementioned biophysical parameters were measured (see Data section for details). Figure 2 shows the joint evolution of HSIC and RMSE. The color intensity along the paths represents the increasing (decreasing) value of the consistency hyperparameter µ as it goes from reducing consistency with the ancillary variable (positive) to increasing consistency (negative). Optimal solutions, i.e. predictions with lower RMSE and higher HSIC to the ancillary variable, are obtained for negative values of µ, meaning that increasing the consistency between the variables (LAI-fCOVER and Chla-LAI) also helps in reducing the RMSE of the prediction. Similarly to the example of parametric models, higher consistency allows for higher accuracy (lower error) up to certain values starting from which a trade-off appears between consistency and accuracy. This example illustrated that, while many accurate solutions are possible when retrieving biophysical parameters with ML, one can achieve an optimal solution that is both accurate and ensures consistent results across variables.

Zoom In Zoom Out Reset image size

Separating forced and internal variability is a key challenge in climate science [34, 38, 39], and in particular in the context of detection and attribution (D&A) of externally forced climate change [40, 41]. D&A typically uses fingerprints that encapsulate the structure of forced warming from model simulations into a spatial or spatiotemporal pattern [40–42]. Observations and climate model control simulations are then projected onto those fingerprints to assess whether a signal, such as multidecadal temperature trends at the global or continental scale, exceed internal climate variations [42].

In recent research, forced climate signals have been extracted from observations via signal-to-noise optimized fingerprints from climate models obtained through statistical learning [39, 43–45], thus building on ideas from traditional D&A. Statistical learning algorithms in this context aim to minimize the influence of internal variability and model disagreement to extract a forced signal [46]. However, patterns of internal variability (such as El Niño Southern Oscillation, ENSO) may nonetheless project onto the fingerprint, for example, if the fingerprint is derived from a given model (or a given set of models) but then applied to an unseen test model (or observations) with a potentially different representation of key modes of internal variability, such as ENSO. This phenomenon may lead to an inaccurate extraction of forced signals from models or observations, and could lead to over- or under-confidence in D&A statements if models systematically under- or overestimate the magnitude of internal variability [44, 47]. Hence, robustness in fingerprint extraction with respect to model structural differences, or potential systematic biases in the models' representation of internal variability or of forced patterns, is an important issue that still remains a key uncertainty. This issue is related to distributional robustness in statistics [48] and transfer learning in ML [49]. For example, climate models simulate sea surface temperature patterns and warming that are broadly consistent with observations on a global scale [50], and also capture key modes of internal variability such as ENSO [51]. However, models differ (1) in the time scales of irregular ENSO fluctuations, (2) in the characteristic patterns of ENSO-related climate variability, and (3) in how ENSO variability is projected to evolve in a warming climate [51]. Moreover, differences across models in the forced warming and variability in multidecadal trends are particularly pronounced in the equatorial Pacific [50, 52]. While historical model simulations show generally warming sea surface temperatures in the equatorial Pacific (e.g. ENSO3.4 region: $[-170^\circ \,\textrm{E}, -120^\circ \,\textrm{E}]$ and $[-5^\circ \,\textrm{N}, 5^\circ \,\textrm{N}]$), observations show only very modest warming in this region over multi-decadal time scales [52]. It is currently unclear whether this discrepancy is due to an unusual realization of the observed climate (i.e. internal variability), or whether the models show systematic biases in the forced response compared to observations [50].

Here, we illustrate the idea of estimating the forced climate response from individual ensemble members (similar to e.g. [44, 46]) but in a way that takes into account the ENSO-related variability by reducing consistency to the ENSO3.4 region variability expressed by the Niño3.4 index. Reducing consistency in this context thus implies that ENSO-related variability does not project onto the fingerprint (i.e. the prediction of the forced response should be independent of the ENSO3.4 region variability), and hence our goal is to attain improved robustness in the context of uncertainties related to variability in the equatorial Pacific.

For this purpose, we employ the multi-model large ensemble archive [34] and the climate variability diagnostics package [35] (more details in Data section). Winter (December–January–February, DJF) seasonal average values of the forced response in each climate model are calculated as the ensemble average of multiple model runs that differ from each other by a small perturbation. By taking the average of the ensemble, the small oscillations between models that are attributed to internal variability are smoothed out and the end product reflects the forced climate response in the absence of variability [35]. These values constitute the target variable $\mathsf{y}$ of our prediction. The statistical model is based on DJF anomalies in individual ensemble members as predictor variables $\mathsf{x}$. We run the PKL with Niño3.4 as the ancillary variable $\mathsf{s}$ that we want the prediction to be independent to (i.e. PKL is run with positive increasing values of the consistency hyperparameter, µ > 0). For illustration, we employ a different climate model for training (CanESM2), validation (CSIRO-Mk3-6-0) and testing (CESM1-CAM5).

Results are presented in figure 3. Without reducing consistency (µ = 0), the prediction $\hat{\mathsf{y}}$ and Niño3.4 are fairly correlated (Pearson's correlation $\rho = -0.72$) and dependent from the prediction for the unseen test model (figure 3(a)). This implies that ENSO-related uncertainties or potential systematic biases may alias into the prediction of the forced response for an unseen model or in observations. However, if Niño3.4 is taken as the ancillary variable and upon reducing consistency, µ > 0, the prediction becomes more decorrelated ($\rho(\hat y,s) = 0.03$) and independent (HSIC $(\hat y,s) = 0.0005$) of the variability introduced by ENSO3.4 (figure 3(b)). The results imply that our forced response estimate would be more robust to ENSO-related uncertainties with an appropriate magnitude of µ across different climate models or between models and observations. An illustration of the effect the independence constraint has for the prediction can be seen in the supplementary material S2.

Zoom In Zoom Out Reset image size

In summary, we have illustrated that estimating a climate signal (such as the forced response) from climate variables can be learnt in a way that is 'consistent' to certain known model structural uncertainties such as the representation of internal variability and warming in the equatorial Pacific that determine ENSO-related variability in the climate system.

The PKL functional for the function class f used in this work, is defined as:

$$\begin{align*} f_\ast & = \mathop{\mathrm{arg\,min}}_{\,f\in\mathcal H} \bigg\{\frac{1}{n}\sum_{i = 1}^n V(\,f(\mathbf{x}_i),y_i) + \frac{\lambda}{n} \Vert \,f\, \Vert _{\mathcal H_k}^2 \nonumber\\ &\qquad\qquad\quad+ \mu \widehat{\text{HSIC}}_{k,l}(\,f(\mathsf{x}),\mathsf{s})\bigg\}, \end{align*}$$

where we adopted the reproducing RKHS $\mathcal H_k$ as the hypothesis class. Given that the selected consistency penalty term, HSIC, only depends on the unknown function f through its evaluations at the training inputs $\{\mathbf{x}_i\}$, direct application of the Representer's theorem [55] tells us that the optimal solution can be written as a linear combination of kernel function evaluations $f = \sum_{i = 1}^n\alpha_i k(\cdot,\mathbf{x}_i)$. Hence, we obtain the so-called dual problem:

$$\begin{align*} \min_{\boldsymbol{\alpha} \in \mathbb{R}^n} \bigg\{\|\mathbf{K}\boldsymbol{\alpha} - \mathbf{y}\|_2^2 + \frac\lambda n \boldsymbol{\alpha}^\top \mathbf{K} \boldsymbol{\alpha} + \frac \mu {n^2} \boldsymbol{\alpha}^\top \mathbf{H}\mathbf{L}\mathbf{H}\mathbf{K} \boldsymbol{\alpha}\bigg\}, \end{align*}$$

which can now be solved for α directly. In the case of the squared L₂ norm for $V: = \|f(\mathbf{x})-\mathbf{y}\|_2^2$, the above dual problem has an analytic solution [27]:

$$\begin{align*} \boldsymbol{\alpha} = \bigg(\mathbf{K} + \frac{\lambda}{n}\mathbf{I} + \frac \mu {n^2} \mathbf{H}\mathbf{L}\mathbf{H}\mathbf{K}\bigg)^{-1} \mathbf{y}. \end{align*}$$

The dependence estimator is sensitive to the hyperparameters of the kernel functions k and l. To avoid variations in HSIC exclusively due to changes in these parameter values, some form of normalization is needed. Nevertheless, one cannot normalize both kernels arbitrarily since the optimization problem would become nonlinear with α and would not admit a closed-form solution anymore. As the parameters from K are adjusted to the data during the optimization process, one could partially normalize the cross-covariance on L. We introduce a modification inspired by the normalized version of HSIC called NOCCO [56]. This modification simply replaces ${\mathbf{H}\mathbf{L}\mathbf{H}}$ with ${\mathbf{H}\mathbf{L}\mathbf{H}}({\mathbf{H}\mathbf{L}\mathbf{H}} + n\epsilon\mathbf{I})^{-1}$, where is a regularization parameter used in the same way as [56]. The solution then becomes:

$$\begin{align*} \boldsymbol{\alpha} = \bigg(\mathbf{K} + \frac{\lambda}{n}\mathbf{I} + \frac \mu {n^2} \mathbf{H}\mathbf{L}\mathbf{H}({\mathbf{H}\mathbf{L}\mathbf{H}} + n\epsilon\mathbf{I})^{-1}\mathbf{K}\bigg)^{-1} \mathbf{y}. \end{align*}$$

In all our experiments, we used the squared exponential (SE) kernel for both k and l, that is e.g. $k(\mathbf{a},\mathbf{b}) = \exp(-\|\mathbf{a}-\mathbf{b}\|^2/(2\sigma^2)),$ which captures sample similarity well in most of the problems, and only one hyperparameter σ needs to be tuned. For the standard kernel ridge regression (KRR) we implemented a standard cross-validation scheme to determine the ridge regression parameters, σ and λ, setting the squared error loss as the metric to minimize. For the ancillary variable we set the SE hyperparameter to the average Euclidean distance between all points, which is a commonly used heuristic in the kernel methods literature [23, 57]. A much larger sigma would make HSIC approximate a linear dependence estimate, while a much smaller legnthscale σ would make HSIC insensitive to dependence [29].

Using HSIC in a kernel-based regression framework can be actually seen as a modified Gaussian process prior. From a Gaussian Process (GP) perspective, the prior covariance corresponds to a posterior covariance for meta-observations encouraging that the predictive function f remains constant as the ancillary variables vary. Let us assume that the loss corresponds to the negative conditional log-likelihood in some probabilistic model, i.e. that $V(\,f(\mathbf{x}_i),y_i) = -\log p\left(y_i|f(\mathbf{x}_i)\right)$, which is true for a wide class of loss functions. We write $V(\mathbf{y},\mathbf{\,f}) = -\sum_{i = 1}^n \log p(y_i|f(\mathbf{x}_i))$ to denote the rescaled conditional negative log-likelihood. Consider now using explicit feature mapping $\mathbf{x}_i\mapsto \phi(\mathbf{x}_i)$ and denoting the feature matrix by Φ, we have $\mathbf{\,f} = \boldsymbol{\Phi}\boldsymbol{\beta}$ and recast optimization as:

$$\begin{align*} \min_{\beta \in \mathbb{R}^m} \bigg\{\frac{1}{\lambda}V(\mathbf{y},\boldsymbol{\Phi}\boldsymbol{\beta}) + \boldsymbol{\beta}^\top \boldsymbol{\beta} + \delta \boldsymbol{\beta}^\top \boldsymbol{\Phi}^\top \mathbf{H}\mathbf{L}\mathbf{H} \boldsymbol{\Phi} \boldsymbol{\beta}\bigg\}, \end{align*}$$

where β play the same role as α in the PKL, and we defined $\delta = \mu/\lambda n$ for convenience. The two regularization terms correspond, up to an additive constant, to a negative log-prior of $\boldsymbol{\beta}\sim\mathcal{N}\left(0,\left(\mathbf{I}+\delta \boldsymbol{\Phi}^\top \mathbf{H}\mathbf{L}\mathbf{H}\boldsymbol{\Phi}\right)^{-1}\right)$, which in turn gives a prior on the evaluations:

$$\begin{align*} \mathbf{\,f}\sim\mathcal{N}\left(0,\boldsymbol{\Phi}\left(\mathbf{I} + \delta \boldsymbol{\Phi}^\top \mathbf{H}\mathbf{L}\mathbf{H}\boldsymbol{\Phi}\right)^{-1}\boldsymbol{\Phi}^\top\right). \end{align*}$$

By directly applying the Woodbury-Morrison formula, the covariance matrix in this prior becomes $(\mathbf{K}^{-1}+\delta \mathbf{H}\mathbf{L}\mathbf{H})^{-1}$, compared to K in the standard GP case. Thus, adding an HSIC regularizer corresponds to modifying the prior on function evaluations f. The posterior mode in a Bayesian model using a modified GP prior becomes:

$$\begin{align*} f \sim \mathcal{GP}\left(0, k(\cdot,\cdot)-k_{\mathbf{X}\cdot}^\top(\mathbf{K}\mathbf{H}\mathbf{L}\mathbf{H}+ \delta^{-1}\mathbf{I})^{-1}\mathbf{H}\mathbf{L}\mathbf{H} k_{\mathbf{X}\cdot}\right). \end{align*}$$

where $k_{\mathbf{X}\cdot} = [k(\cdot,\mathbf{x}_1),\cdots,k(\cdot,\mathbf{x}_n)]^\top$, for any training set $\{\mathbf{x}_i\}_{i = 1}^n$. Treating the learning problem in the GP framework actually allows to derive uncertainty estimates and hyperparameter tuning using marginal likelihood maximization.

Neural networks can also benefit from the new loss including the dependence term. Including the HSIC term in standard backpropagation leads to simple updating rules for training neural networks. Let us denote a feedforward neural network function $\hat{\mathbf{y}} = {\mathcal F}(\mathbf{X};\mathbf{W})$, parameterized by a set of weights W. The backpropagation algorithm uses gradient descent to adjust model weights iteratively $\mathbf{W}[t]\leftarrow\mathbf{W}[t-1] + \eta\nabla J_W$, where $J = \|\mathbf{e}\|$ is the loss (cost, energy) function that depends on training error $\mathbf{e} = \mathbf{y}-\hat{\mathbf{y}}$ and η > 0 is the learning rate. Including the physics-dependence loss in a neural network training is straightforward, and only involves replacing the output error $\mathbf{e} = \mathbf{y} - \hat{\mathbf{y}}$ with $\mathbf{e}^{^{\prime}} = \mathbf{y} - (\lambda\mathbf{I}+\mu\tilde{\mathbf{S}} \tilde{\mathbf{S}}^\top)\hat{\mathbf{y}}$, which intuitively trades the training error for physics consistency. The standard backpropagation algorithm can be applied to the new error function, $J^{^{\prime}}_W: = \|\mathbf{e}^{^{\prime}}\|^2$. Alternatively one could use other training algorithms like Adam or automatic differentiation if the interest is to learn the dependence measure parameters end-to-end. The algorithm allows us to train large scale problems while encoding physical consistency concepts easily.