Can Neural Networks Think Like Geomorphologists?
This work was recently published in Geophysical Research Letters. This page contains a more casual write-up of the key parts of the work. The paper (and supplement) contains a more complete write-up.
One of the motivating questions for my work is how we can quantify important but complex patterns in spatial data. Even relatively complex metrics can’t capture all of the relevant information in spatial data. Neural networks are a promising tool to help find and quantify complex spatial patterns, but as an earth scientist, I want to be able to tie machine learning results back to physical theory in a way that enhances our understanding of the natural world.
So as a proof of concept, here’s some work where we trained a simple neural network to invert topography generated by a simple landscape evolution model, and then dug into the model structure to understand what geomorphology it has (and has not) learned.
The streampower-diffusion model
For our simple generated landscapes we use the streampower-diffusion model which generates landscapes with the parameter \(D/K\) as a key control. The model domain and parameters were inspired by Perron et al. 2008.
Streampower-diffusion model details (click to expand)
We chose to use landscape evolution models in this project, because we wanted to bypass the issues and complexities with real data. By using this synthetic data we get to focus on evaluating a neural network’s learning of a simple, highly controlled, and well understood system. Particularly we use the “streampower-diffusion” (SPD) model, which is a very simple and very commonly used model for landscape evolution. The general idea is that in (some) landscapes there are two primary sources of erosion: overland flow and hillslope processes. Erosion due to overland flow is advective, and erosion due to hillslope processes is diffusive. Geomorphologists have developed some physics-based models that quantify this, which we won’t go into, but that give us this equation for erosion: \(\frac{\partial z}{\partial t}=D\nabla^2z - KA^m|\nabla z|^n+U\) The key things you should know are that \(\frac{\partial z}{\partial t}\) really means the rate of change of elevation (so basically erosion). The first part is hillslope diffusion, and \(D\) controls how effective diffusion on the landscape is. The second part is fluvial erosion (the streampower law) and \(K\) controls how effective that is. Lastly, \(U\) is uplift due to tectonics (assumed to be constant in this system). In this project we vary \(D\) and \(K\), solve this equation numerically to simulate 3 million years of landscape evolution, and then play with the results.
I’ve mentioned advection and diffusion, and so if you have a certain type of background, your mind might be screaming “Péclet number”, and yes, the ratio of advective to diffusive transport is a key control on this system. In this work we look at the ratio \(D/K\) which can be thought of as a sort of inverse Péclet number. The reason for \(D/K\) is a since-aborted attempt to compare with a different data set, but we initially did this work looking at \(K/D\) (this comes back in the paper where we notice an interesting diversion between input preference for \(D/K\) and \(K/D\)) and non-dimensionalized Péclet numbers, and it all works about the same. By keeping our ratio dimensionalized we avoid having to choose length scales which ends up requiring some justification that really isn’t necessary for the core results. I discuss this a little bit more here.
So let’s look at how changing \(D/K\) changes the landscape that this model produces:
Ok so it’s pretty clear that landscapes with different \(D/K\) values look different in a systematic way. And I bet that you would be able to distinguish between a high and low \(D/K\) landscape, but can you systematically calculate \(D/K\) from bulk topography? We know that \(D/K\) controls the spacing of valleys Perron et al. 2008, and the drainage density Howard 1997Tucker and Bras 1998, but these relationships aren’t necessarily one-to-one. A convolutional neural network can predict numbers like \(D/K\), and does so by extracting spatial features from the inputs.
The neural network
Some neural network basics (click to expand)
The neural network is a function that takes in a 2D single-band input, and produces a single number. It does so by feeding the input through a series of “layers”, or simple functions. The network has four types of layers: convolutions, rectified linear units (ReLU), maximum pooling layers (MaxPool) and linear layers.
- Convolution
- a convolutional layer contains a window of weights (kernel) which “sweeps” over the input. For each section of input the window overlaps with, the convolution outputs the weighted sum of the input, with each pixel being modified by the weight overlaying it. The convolutional layers are the parts of the network that are most sensitive to spatial features. Traditionally the convolution is thought to be sensitive to a specific feature, and the output is referred to as a “feature map”, with the value of each pixel of output representing the level of detection of the feature the convolution is sensitive to at that point of the input.
- ReLU
- a ReLU (rectified linear unit) is a simple ramp function that sets the output to zero if below zero, and leaves it unchanged if zero or above. The role of this function is to introduce nonlinearities into the network.
- MaxPool
- This is functionally a resampling layer that takes an input, divides it into subsets (in this network, cells of two-by-two pixels), and returns only the maximum of that subset. This helps reduce the size of the data and network.
- Linear
- These are \(n\) by \(m\) linear transformations that help reduce the final feature maps into a single number that is interpreted as the \(D/K\) ratio.
Each convolution has a number of channels, which is the number of kernels that transform the input. In the network used in this paper, the first convolutional layer has ten channels, meaning there are ten three-by-three kernels that transform the input into ten different feature maps, or a ten-band image. Each output feature map is then subjected to the ReLU operation, and the maxpooling operation. This final ten-band output is the input to the next set of convolutions. The second convolutional layer has five-by-five convolutions and 20 channels, meaning that there are 20 distinct kernels, each of which is applied to the entire ten-band input. This means each kernel in this layer is a five-by-five-by-ten window of weights. Again, these 20 outputs are subjected to ReLU and maxpooling. The final convolutional/ReLU/MaxPool layer is a 20-band seven-by-seven convolution. This 20-band 2D output is flattened into a (very long) 1D vector. This is what is passed to a series of linear transformations and ReLU layers to become a single number. Specifically, it is passed through a linear transformation with 100 outputs and subjected to a ReLU activation layer. This is then passed to a linear transformation and ReLU with ten outputs, and finally a linear transformation that produces a single number output that is interpreted as the network’s inferred \(D/K\) value for the given input.
We use a pretty basic architecture with three convolution layers followed by three linear layers. See this page for a diagram of the neural network (it includes feature maps and activation maximization images).
The question
So here’s the question: Can a neural network extract \(D/K\) from bulk topography? And if so, does it use geomorphically meaningful features to do so? Let’s run thousands of models and see if a neural network can invert them.
The answer
It’s yes (to the first part).
There is some serious spread at low \(D/K\) values but there’s a lot less data there just due to how we chose parameter spacing. Also look above at how the landscapes change as \(D/K\) is lowered, they get more and more dissected. It’s possible that landscapes are really really similar in that low \(D/K\) range, where the scale of dissection is smaller than the scale of our model resolution.
Still, we have a normalized root mean squared error of 0.0207 which I think is good enough to say a neural network can extract \(D/K\) from bulk topography. And it does about an order of magnitude better when trained on slope and curvature. So is it doing so in a geomorphically meaningful way? Is it thinking like a geomorphologist?
We tried getting at this in a lot of different ways (see this blog post) but there are two tests that were the most successful. I said before that we know that \(D/K\) controls the spacing of valleys Perron et al. 2008, and the drainage density Howard 1997Tucker and Bras 1998, and we also know the general shape of these relationships: more streampower should mean narrower valleys and higher drainage density, so let’s test if the neural network has this same understanding. We do this by creating two datasets: an “all valley spacing” and an “all drainage density” dataset. The “all valley spacing” dataset is a corrugated sine surface where elevation is just \(\sin(y)\) for some periodicity. The resulting landscape is really just a bunch of parallel “valleys” with exact spacing. So when we take the neural network trained to infer \(D/K\) from the streampower-diffusion model topography and give it these sine landscapes, what does it find? Well, that a landscape with larger valley spacing should have a larger \(D/K\) value.
We did a similar test with drainage density. We created a bunch of drainage network masks (value 0 where there is a channel, value 1 where there is none), functionally a flat plane with incised channels, nothing else.
Here we see drainage density decrease as \(D/K\) increases, meaning that there is relatively less stream erosion on the landscape. We talk a bit more about the weird breaks in the relationship in the paper. These are both generally (as in the right direction) the relationships we as geomorphologists would expect! This is, in my opinion, the most interesting result of the paper and the best evidence that the neural network has a “geomorphic understanding”.
What the network has failed to learn
The above result is exciting, but it doesn’t mean this neural network is some generalized \(D/K\) machine. If you look at the model runs above, you’ll probably notice that they’re all really similar in a way that real landscapes are not, and it will probably not surprise you to hear that this network can’t even discriminate between two distinct real landscapes when given real data. We did some tests to try to get a better sense of how and where it fails. One way we did this is by taking the network that was trained on “normal” landscape evolution data, and seeing how it performed when given inputs that were swapped or shuffled around. We split the image in half and swapped it to make a central ridge instead of a central valley, and broke the image into nine tiles and shuffled those around. There are some example inputs in the paper supplement.
The network trained on elevation data totally fails to generalize, meaning that it really needs the global structure we imposed in our model run domain to be able to perform. Interestingly, slope does not have that same issue. So not only do geomorphic derivatives allow a neural network to perform better on this task than the elevation data, but they also help it learn the task in a more robust way! This perhaps isn’t so surprising since slope and curvature show up in the governing equation of the model, and are more locally distributed than elevation. We discuss this more in the paper and the supplement which I urge you to check out.
What’s next?
Well, there are some rabbit holes about the way that performance is impacted by the choice of data type and target that should be chased down, and some natural extensions to more complex model domains, more complex geomorphic problems, and real data. I hope to have more on all of those soon!