physics-informed invertible neural networks for hydrodynamic imaging

Note The work presented here is an attempt at making a difficult topic easily understandable and does not represent proper scientific rigour in every aspect. For scientific rigour, I would like to refer you to the original thesis.

A project where I spent the better part of a year researching fish, physics, and neural network architectures. The setup is a small sphere moving through water while an array of eight sensors on the bottom records flow disturbances created by the moving object.

This scene captures the essence of hydrodynamic imaging, which is the use of water-motion patterns to detect and characterize nearby objects. Fish developed this ability over millions of years through their lateral line organs, and current research tries to replicate that capability in machines.

The experimental setup is simple. Eight optical sensors are arranged in a line across the bottom of a virtual pool, each 50 millimeters apart, creating a 400-millimeter sensing array. Above them, a sphere follows several trajectories, including straight, curved, and zigzag paths.

As the sphere moves, it creates predictable disturbances in the water that follow potential flow theory. Each movement generates a distinct flow signature that the sensors can detect.

Understanding Potential Flow: The Invisible Architecture of Moving Water

Before diving into the mathematics, it helps to clarify what potential flow represents. A useful mental model is a three-dimensional grid of water particles. When the sphere moves through this grid, nearby particles adjust their motion according to physical constraints.

This adjustment follows strict rules. Water is treated as incompressible, the flow is irrotational, and viscosity is neglected. These constraints create structured patterns. Picture dropping a stone into a still pond. Ripples spread outward in circles because the water follows mathematical laws.

In three dimensions around the moving sphere, streamlines wrap around the object, squeeze together near its widest section, and spread apart elsewhere. These streamlines carry information: density indicates flow speed, direction indicates motion, and the overall pattern relates to object size.

⚡ Key Insight: Fish infer nearby objects from changes in local flow patterns measured by the lateral line.

When the sphere moves through the pool, it creates patterns governed by mathematical laws. A central concept is the velocity potential, a scalar function φ(x,y,z) that describes the flow around the sphere.

For incompressible, irrotational flow (which is a reasonable approximation for our sensing scenario), we can define the velocity field as the negative gradient of this potential:

$$\vec{v} = -\nabla \phi = -[\partial\phi/\partial x, \partial\phi/\partial y, \partial\phi/\partial z]$$

This means if we know φ everywhere, we know the entire flow field. For a sphere of radius a moving with velocity $\vec{W}$ through an inviscid fluid, the velocity potential becomes:

$$\phi(x,y,z) = (a^3/2) \times \frac{(\vec{W} \cdot \vec{r})}{ \lvert\vec{r}\rvert^3}$$

where $\vec{r}$ = [x,y,z] is the position vector from the sphere’s center. This equation shows that the flow disturbance depends on the cube of sphere size ($a^3$), scales linearly with velocity ($\vec{W}$), and falls off with the cube of distance ($\vert\vec{r}\vert^3$).

An important consequence is that this velocity potential must satisfy Laplace’s equation:

$$\nabla^2\phi = \frac{\partial^2\phi}{\partial x^2} + \frac{\partial^2\phi}{\partial y^2} + \frac{\partial^2\phi}{\partial z^2} = 0$$

This is a fundamental physical constraint that any realistic flow must satisfy. The flows we detect are solutions to this partial differential equation, which means they carry the fingerprints of the underlying physics.

For our sensor array positioned along the x-axis, we can derive the specific velocity components that our sensors detect. After some calculus, the flow patterns decompose into three fundamental “wavelets”:

def compute_wavelets(rho):
    """
    Compute the three fundamental wavelets for lateral line sensing

    Args:
        rho: normalized sensor position (s-x)/y where s is sensor position,
             (x,y) is sphere position

    Returns:
        even_wavelet, odd_wavelet, navelet
    """
    denominator = (1 + rho**2)**(5/2)

    even_wavelet = (1 - 2*rho**2) / denominator      # Ψₑ: parallel motion
    odd_wavelet = -3*rho / denominator               # Ψₒ: perpendicular motion
    navelet = (2 - rho**2) / denominator             # Ψₙ: directional info

    return even_wavelet, odd_wavelet, navelet

These wavelets are the mathematical building blocks of hydrodynamic sensing. The even wavelet captures signatures of parallel motion, the odd wavelet captures perpendicular motion, and the navelet provides additional directional information.

The velocity components our sensors actually measure combine these wavelets in specific ways:

$$v_x = \frac{a^3 \left\lVert w \right\rVert}{2d^3}(\Psi_o \sin{\theta} - \Psi_e \cos{\theta}),$$ $$v_y = \frac{a^3 \left\lVert w \right\rVert}{2d^3}(\Psi_n \sin{\theta} + \Psi_o \cos{\theta}).$$

Here, $d$ is the distance from sphere to sensor array, and $\theta$ is the angle of motion. Sphere size appears as $a^3$, so a sphere twice as large creates eight times the signal. Distance also appears cubed in the denominator, so doubling the distance reduces the signal by a factor of eight. This cubic relationship makes detection of small, distant objects challenging.

The image shows a graphical representation of the wavelet family (a) and the velocity profiles, vx (b) and vy (c) for a sphere of size 10 (mm) moving parallel to the array at a constant speed of 10 mm/s. This is the concatenated representation of 1024 uniformly spaced sensors between -200 mm and 200 mm.

def compute_velocity_profiles(sphere_pos, sphere_velocity, sphere_radius,
                            sensor_positions, angle):
    """
    Compute velocity profiles at sensor locations for a moving sphere
    """
    x, y = sphere_pos
    sensors = np.array(sensor_positions)

    # Compute normalized positions for each sensor
    rho = (sensors - x) / y
    distance = np.sqrt(x**2 + y**2)

    # Get wavelets for all sensor positions
    psi_e, psi_o, psi_n = compute_wavelets(rho)

    # Velocity magnitude scaling factor
    velocity_magnitude = np.linalg.norm(sphere_velocity)
    scale_factor = (sphere_radius**3 * velocity_magnitude) / (2 * distance**3)

    # Compute velocity components
    v_x = scale_factor * (psi_o * np.sin(angle) - psi_e * np.cos(angle))
    v_y = scale_factor * (psi_n * np.sin(angle) + psi_o * np.cos(angle))

    return v_x, v_y

These mathematical relationships connect moving spheres to sensor measurements through well-understood equations. Fish use similar flow cues through the lateral line. This project applies those same cues in a machine-learning setting.

So what can we extract from these hydrodynamic signatures? Traditional research in this field focused on two outputs: object location (x, y) and heading angle ($\theta$). The equations also contain information about sphere size ($a^3$) and speed ($\lVert\vec{W}\rVert$). A key question was whether those additional characteristics could be recovered from sensor data.

This became one of the central challenges of my research. Previous work had established that location and angle were recoverable, while speed and size had not been extracted successfully from artificial lateral line data. The mathematical relationships were clear, and the inverse problem, which works backward from noisy sensor measurements to object properties, remained difficult.

Working backwards from sensor readings to object properties is what mathematicians call an inverse problem. It is often ill-posed, meaning multiple solutions can exist or small measurement changes can produce large output changes. In our case, consider the equation:

$$signal\ strength \propto (a^3 \times \lVert\vec{W}\rVert) / d^3$$

A small, fast sphere close to the sensors (small $a$, large $\lVert\vec{W}\rVert$, small $d$) might create exactly the same signal as a large, slow sphere further away (large a, small $\lVert\vec{W}\rVert$, large $d$). This is the fundamental ambiguity that makes inverse problems so challenging.

⚠️ The Fundamental Challenge: Inverse problems are like trying to figure out what someone looks like from their shadow. Multiple different shapes can cast very similar shadows, especially when the shadow is blurry or distorted.

Real sensors add another layer of complexity: noise. The measured velocities include random error:

$$v_\text{measured} = v_\text{true} + n\ \text{where}\ n \sim N(0, σ²)$$

The signal-to-noise ratio (SNR) becomes:

$$SNR = \frac{\lvert v_\text{true}\rvert}{\sigma} = \frac{(a^3 \times \lVert\vec{W}\rVert \times \text{wavelet amplitude})}{ (d^3 \times \sigma)}$$

This means small, slow, or distant objects become exponentially harder to detect and characterize. When SNR drops below about 3, reliable detection becomes nearly impossible.

Another approach is to train neural networks to learn the mapping from sensor data to object properties directly from examples. Previous research had established that Long Short-Term Memory (LSTM) networks perform well on this task because temporal context carries important information.

While a single snapshot of sensor data might be ambiguous, watching how the signals evolve over time provides crucial additional information. LSTMs are a specialized type of recurrent neural network designed to capture long-term dependencies in sequential data through sophisticated gating mechanisms that control information flow and memory retention. For a detailed explanation of LSTM architecture and functionality, see the Wikipedia article on LSTMs.

class LateralLineLSTM(tf.keras.Model):
    def __init__(self, hidden_units=256, window_size=16):
        super().__init__()
        self.window_size = window_size
        self.lstm = tf.keras.layers.LSTM(
            hidden_units,
            activation='tanh',
            return_sequences=False
        )
        self.dense = tf.keras.layers.Dense(3)  # x, y, theta

    def call(self, inputs):
        # inputs shape: (batch_size, time_steps, features)
        # features = 16 (8 sensors × 2 velocity components)
        lstm_out = self.lstm(inputs)
        predictions = self.dense(lstm_out)
        return predictions

A single measurement can be ambiguous, while a short sequence provides enough context to identify the pattern. The LSTM uses that temporal context to interpret hydrodynamic signatures.

I also evaluated more sophisticated architectures. The first was Invertible Neural Networks (INNs), a class of networks that can learn bidirectional mappings. Traditional models usually learn X → Y. INNs learn both X → Y and Y → X.

This works through “coupling layers” that implement a specific mathematical structure. For an input vector x split into two parts x₁ and x₂:

def coupling_layer_forward(x1, x2, scale_net, translation_net):
    """
    Forward pass through a coupling layer
    """
    s = scale_net(x1)      # Scale parameters
    t = translation_net(x1)  # Translation parameters

    y1 = x1                           # Pass through unchanged
    y2 = x2 * tf.exp(s) + t          # Affine transformation

    return y1, y2

def coupling_layer_inverse(y1, y2, scale_net, translation_net):
    """
    Inverse pass through the same coupling layer
    """
    s = scale_net(y1)      # Same scale network
    t = translation_net(y1)  # Same translation network

    x1 = y1                           # Direct recovery
    x2 = (y2 - t) * tf.exp(-s)       # Exact inverse

    return x1, x2

🔄 Invertibility: Half of the input passes through unchanged while the other half undergoes a designed transformation. Because the transformation is known, the inverse can be computed exactly.

Half of the input passes through unchanged, while the other half undergoes an affine transformation (scaling and translation) based on the unchanged half. Because affine transformations are easily invertible, we can always recover the original input exactly.

The Jacobian determinant (needed for proper probability calculations) has a simple form:

$$\lvert \text{det}(\frac{\partial y}{\partial x}) \rvert = \exp(\sum_i s_i)$$

where sᵢ are the outputs of the scale network. This exponential form is both easy to compute and numerically stable.

For our lateral line problem, we can train the network to learn the mapping from object properties (x, y, $\theta$) to sensor readings ($v_x$, $v_y$), and we automatically get the inverse mapping for free:

class LateralLineINN(tf.keras.Model):
    def __init__(self, num_coupling_layers=4, hidden_dims=64):
        super().__init__()
        self.coupling_layers = []

        for i in range(num_coupling_layers):
            # Each coupling layer has its own scale and translation networks
            scale_net = tf.keras.Sequential([
                tf.keras.layers.Dense(hidden_dims, activation='relu'),
                tf.keras.layers.Dense(hidden_dims, activation='relu'),
                tf.keras.layers.Dense(hidden_dims, activation='relu'),
                tf.keras.layers.Dense(hidden_dims, activation='relu'),
                tf.keras.layers.Dense(16)  # Output dim for half the features
            ])

            translation_net = tf.keras.Sequential([
                tf.keras.layers.Dense(hidden_dims, activation='relu'),
                tf.keras.layers.Dense(hidden_dims, activation='relu'),
                tf.keras.layers.Dense(hidden_dims, activation='relu'),
                tf.keras.layers.Dense(hidden_dims, activation='relu'),
                tf.keras.layers.Dense(16)  # Output dim for half the features
            ])

            self.coupling_layers.append((scale_net, translation_net))

    def forward(self, x):
        # x: object properties [x, y, theta] + zero padding
        for scale_net, translation_net in self.coupling_layers:
            x = self.coupling_step_forward(x, scale_net, translation_net)
            x = self.permute(x)  # Permute dimensions between layers
        return x  # sensor readings [v_x, v_y] + latent variables

    def inverse(self, y):
        # y: sensor readings + latent variables
        for scale_net, translation_net in reversed(self.coupling_layers):
            y = self.permute_inverse(y)
            y = self.coupling_step_inverse(y, scale_net, translation_net)
        return y  # object properties

There is a subtlety here. Real inverse problems are often ill-posed, meaning multiple inputs can produce the same output. INNs handle this by augmenting the output space with latent variables $z$ sampled from a known distribution:

$$\text{Forward}: (x, y, \theta) + \text{zero padding} \rightarrow (v_x, v_y) + \text{latent variables}$$ $$\text{Inverse}: (v_x, v_y) + \text{latent variables} \rightarrow (x, y, \theta) + \text{recovered padding}$$

The latent variables capture the “missing information” needed to make the inverse mapping unique. During training, we ensure these latent variables follow a standard Gaussian distribution using Maximum Mean Discrepancy (MMD) loss.

The second innovation drew from Physics-Informed Neural Networks (PINNs). The core idea is to incorporate known physical laws directly into the learning process.

Remember our Laplace equation constraint? For potential flow around a sphere, this translates to a specific relationship between the measured velocities. Starting from the velocity potential and taking the required derivatives, we get:

$$\frac{\partial v_x}{\partial x} + \frac{\partial v_y}{\partial y} + \frac{3(W_xx + W_yy)}{(x^2 + y^2)^{(5/2)}} = 0$$

The image shows a Physics-informed neural network where on the left side, we make use of a “regular” neural network. On the right side we then enhance it with the error that derives from a partial differential equation. This physical law becomes an additional “loss term” during training:

def physics_loss(x, y, v_x, v_y, W_x, W_y):
    """
    Compute the physics-informed loss based on Laplace equation
    """
    # Use automatic differentiation to compute partial derivatives
    with tf.GradientTape() as tape:
        tape.watch([x, y])
        # Predict velocities from current position
        v_x_pred, v_y_pred = model(tf.stack([x, y], axis=1))

    # Compute gradients
    dvx_dx = tape.gradient(v_x_pred, x)
    dvy_dy = tape.gradient(v_y_pred, y)

    # Physics constraint from Laplace equation
    distance_squared = x**2 + y**2
    physics_term = 3 * (W_x * x + W_y * y) / (distance_squared**(5/2))

    # Residual should be zero for physically consistent flows
    residual = dvx_dx + dvy_dy + physics_term

    return tf.reduce_mean(tf.square(residual))

🧠 Physics as Teacher: The network learns from examples and from the governing physical laws.

The network is trained to fit data while respecting fluid-dynamics constraints.

I combined both innovations into Physics-Informed Invertible Neural Networks (PI/INNs). Traditional PINNs enforce physics on the forward direction. With an INN, training the forward direction with physics constraints also provides an inverse mapping that inherits those constraints.

The complete loss function balances multiple objectives:

def total_loss(x_true, y_true, theta_true, v_x_true, v_y_true):
    """
    Combined loss function for PI/INN
    """
    # Forward pass: (x,y,θ) → (v_x,v_y,z)
    forward_output = model.forward(tf.stack([x_true, y_true, theta_true], axis=1))
    v_x_pred, v_y_pred, z_pred = tf.split(forward_output, [8, 8, 32], axis=1)

    # Backward pass: (v_x,v_y,z) → (x,y,θ)
    backward_input = tf.concat([v_x_true, v_y_true, z_pred], axis=1)
    backward_output = model.inverse(backward_input)
    x_recon, y_recon, theta_recon = tf.split(backward_output, [1, 1, 1], axis=1)

    # Data fitting losses
    forward_loss = tf.reduce_mean(tf.square(v_x_pred - v_x_true) +
                                 tf.square(v_y_pred - v_y_true))
    backward_loss = tf.reduce_mean(tf.square(x_recon - x_true) +
                                  tf.square(y_recon - y_true) +
                                  tf.square(theta_recon - theta_true))

    # Physics constraint loss
    physics_loss_val = physics_loss(x_true, y_true, v_x_pred, v_y_pred,
                                   W_x_true, W_y_true)

    # Latent distribution loss (MMD to enforce Gaussian distribution)
    z_target = tf.random.normal(tf.shape(z_pred))
    latent_loss = maximum_mean_discrepancy(z_pred, z_target)

    return forward_loss + backward_loss + 0.1 * physics_loss_val + 0.01 * latent_loss

This creates a network that simultaneously learns accurate mappings, respects physical laws, and handles the ill-posed nature of inverse problems.

With the different neural architectures ready, it was time to put them to the test. I created a comprehensive dataset with spheres of different sizes (10-50mm) moving at different speeds (10-50mm/s) along increasingly complex trajectories. The trajectories ranged from simple straight-line motion to complex sinusoidal paths, creating a dataset of around 80 million data points.

First, I tested all three architectures (LSTM, INN, PI/INN) on the traditional tasks that previous research had tackled. The results told a clear story: the LSTM remained the champion for location estimation, achieving errors of 9-18mm on simple paths and maintaining reasonable performance even on complex trajectories. The PI/INN showed promise, consistently outperforming the plain INN and demonstrating that physics constraints do help, especially on complex paths where the additional structure provides valuable guidance.

My research broke new ground when I tried extracting speed and volume characteristics. Using the primary estimates from each network, I then used these estimates to derive the characteristics analytically:

def extract_speed(positions, timestamps):
    """
    Extract object speed from position estimates over time
    """
    # Compute displacement vectors
    displacements = positions[1:] - positions[:-1]
    time_deltas = timestamps[1:] - timestamps[:-1]

    # Instantaneous speeds
    speeds = np.linalg.norm(displacements, axis=1) / time_deltas

    # Average speed over trajectory (robust to noise)
    return np.mean(speeds)

def extract_volume(v_x, v_y, position, angle, speed, sensor_positions):
    """
    Extract sphere radius using inverse velocity potential
    """
    x, y = position
    best_estimates = []

    for i, sensor_pos in enumerate(sensor_positions):
        # Skip sensors with weak signals (low SNR)
        if abs(v_x[i]) < 1e-6:
            continue

        # Compute normalized position for this sensor
        rho = (sensor_pos - x) / y
        distance = np.sqrt(x**2 + y**2)

        # Get wavelet values
        psi_e, psi_o, psi_n = compute_wavelets(rho)

        # Expected velocity component based on wavelets
        expected_v_x = (psi_o * np.sin(angle) - psi_e * np.cos(angle))

        if abs(expected_v_x) > 1e-6:  # Avoid division by zero
            # Solve for radius: v_x = (a³||W||)/(2d³) × wavelet_combo
            radius_cubed = (2 * distance**3 * v_x[i]) / (speed * expected_v_x)
            radius = np.cbrt(abs(radius_cubed))
            best_estimates.append(radius)

    return np.median(best_estimates) if best_estimates else 0

The mathematics here reveals something about error propagation. Since derived quantities depend on primary estimates, their accuracy is fundamentally limited by the underlying predictions. If our position estimates have error σₓ, our speed estimates will have error approximately:

$$\sigma_\text{speed} \approx \sqrt{2} \times \frac{\sigma_x}{\Delta t}$$

This means speed estimation is inherently more sensitive to noise than position estimation, explaining why good primary estimation is crucial for reliable secondary features.

The results showed that the LSTM achieved speed errors of 0.5-0.7 mm/s and volume errors of 0.2-1.9 mm, while the PI/INN showed moderate performance. The study also demonstrated extraction of object speed and volume from artificial lateral line sensor data, which had not been reported previously in this context.

Real-world sensors are noisy, so I tested how each architecture handled increasingly challenging conditions. Remarkably, in high-noise conditions, the physics-informed approach sometimes performed worse than the plain invertible network.

The reason illuminates an important principle in physics-informed learning. When noise is significant, physics constraints create competing objectives in the loss function:

$$\text{Total Loss} = \text{Data Fitting Loss} + \lambda \times \text{Physics Loss}$$

The data fitting term tries to match noisy observations, while the physics term enforces clean mathematical relationships. When noise is large, these objectives conflict because the noisy data no longer satisfies the physical constraints. The optimization process then balances two competing targets, which can lead to weaker solutions.

⚖️ The Physics-Noise Trade-off: With cleaner data, physics constraints improve training. With noisier data, those constraints can reduce fit quality.

This can be seen mathematically. If our true signal follows the physics constraint $f(x,y) = 0$, but our noisy measurements satisfy $f(x,y) = \epsilon$ where $\epsilon$ is the noise contribution, then enforcing $f(x,y) = 0$ actually pulls us away from the best fit to the noisy data.

The LSTM’s consistent superiority reinforced the critical importance of temporal context in hydrodynamic sensing. While individual velocity profiles theoretically contain complete information, temporal sequences provide noise averaging, consistency checking, and predictive context that dramatically improve accuracy. The temporal averaging effect follows:

$$\text{Effective SNR} \approx \sqrt{N} \times \text{Single Measurement SNR}$$

where N is the number of temporal observations. This square-root improvement with sample size is a key insight from statistics, and it explains why temporal models consistently outperform instantaneous approaches.

Physics-informed constraints help when conditions are favorable, such as clean data, well-understood physics, and interpolation tasks. They can hurt when data is noisy or model assumptions are imperfect. Physics acts as a prior, and its influence must be tuned to match data quality.

Successfully extracting speed and volume revealed the importance of error propagation in derived quantities. The mathematical relationships show that errors in primary estimates (position, angle) propagate through the analytical derivations to affect secondary features. For the velocity potential inversion:

$$\partial(\text{radius})/\partial(\text{velocity}) \propto \frac{1}{\text{wavelet amplitude}}$$

This means errors are amplified when the wavelet amplitude is small, which happens when the object is far from sensors or moving in unfavorable orientations. Understanding these mathematical sensitivities helps explain why some configurations are much easier to characterize than others.

This research sits at the intersection of bio-inspired engineering, physics-aware machine learning, and inverse-problem solving from observations to causes. The extraction of additional features shows that careful problem formulation can extend what existing sensor setups can recover.

Fish still outperform current artificial systems in underwater flow sensing through their lateral lines. This research is one step toward that level of capability by teaching models to interpret subtle hydrodynamic signatures. For more on this topic, see work by Ben Wolf and Sietse van Netten, including this paper.

This work also highlights how biological sensing systems can guide technical design. Natural systems in rivers, oceans, and streams remain a useful source of ideas for future sensing methods.