SkyBathy Theory: Data-Driven Lyzenga Model

The SkyBathy tool implements a physics-based approach, solved statistically, to derive bathymetry (water depth) from multispectral satellite imagery. This method is a generalized implementation of the classic Lyzenga algorithm (1978, 2006).

Unlike pure analytical methods, which require precise field measurements of water optical properties (attenuation coefficients $K_d$) and bottom reflectance ($R_b$), this method estimates these parameters directly from the data by using a regression between pixel values and known depth points (from the uploaded CSV).

1. The Physical Basis

The theory is grounded in the Beer-Lambert Law, which describes the exponential attenuation of light as it travels through water. The radiance measured by the satellite ($L_{sat}$) in a visible band can be simplified as:

$$ L_{sat} = L_{\infty} + (L_{bottom} - L_{\infty}) \cdot e^{-2Kz} $$

Where:

$L_{sat}$: Radiance measured at the sensor.
$L_{\infty}$: Deep water radiance (includes atmospheric scattering and surface reflection, but no signal from the bottom).
$L_{bottom}$: Radiance reflected by the seabed (before attenuation).
$K$: Diffuse attenuation coefficient of the water.
$z$: Water depth.
The factor $2$ accounts for the two-way path of light (down to the bottom and back up to the sensor).

Figure 1: Bathymetric map from multispectral images.

2. Linearization (Lyzenga Transform)

To solve the equation for depth $z$ using linear statistical techniques, we must remove the exponential component.

Step A: Deep Water Correction

We subtract $L_{\infty}$ from both sides to isolate the signal coming from the water column and the bottom. In the code, $L_{\infty}$ is automatically calculated as the average of the darkest pixels (5th percentile) within the water mask.

$$ L_{sat} - L_{\infty} = (L_{bottom} - L_{\infty}) \cdot e^{-2Kz} $$

Step B: Logarithmic Transformation

By applying the natural logarithm ($\ln$) to both sides, the equation becomes linear with respect to $z$:

$$ \ln(L_{sat} - L_{\infty}) = \ln(L_{bottom} - L_{\infty}) - 2Kz $$

Rearranging to isolate $z$:

$$ z = \frac{\ln(L_{bottom} - L_{\infty})}{2K} - \frac{1}{2K} \cdot \ln(L_{sat} - L_{\infty}) $$

3. The "Data-Driven" Approach (Multi-Band)

In a pure theoretical model, we would need to know the physical values of $K$ and the bottom reflectance. However, these vary spatially and are difficult to measure without expensive equipment. The approach implemented in SkyBathy bypasses this problem by assuming the equation above can be written as a Multiple Linear Regression Model.

Instead of solving for a single band, we use $N$ visible bands (Blue, Green, Red) to compensate for variations in bottom type (e.g., sand vs. seagrass). The final equation becomes:

$$ z = a_0 + \sum_{i=1}^{N} a_i \cdot X_i $$

Where:

$z$: Estimated Depth.
$X_i$: The transformed variable for band $i$, calculated as:
$$ X_i = \ln(L_{band\_i} - L_{\infty, i}) $$
$a_0, a_1, \dots, a_n$: Empirical calibration coefficients.

What do the coefficients $a_i$ represent?

Unlike the physical method, we do not fix $K_d$ or $R_b$. Instead:

The coefficients $a_i$ (slopes) implicitly absorb the term $-\frac{1}{2K}$ for each band.
The coefficient $a_0$ (intercept) absorbs the constants related to the average bottom reflectance and residual atmospheric effects.

These coefficients are determined by solving a system of linear equations (Ordinary Least Squares) using the control points provided in your CSV file.

4. Implementation in Code

In the file skybathy_derived.html, the process occurs in the solveLyzenga() function following these steps:

Estimate $L_{\infty}$: The function computeDeepWaterRadiance() finds pixels where water is deep (minimum signal) and calculates the subtraction value for each band.
Build Matrix X: For each point in the CSV (x,y,z), the code extracts the pixel value in Blue, Green, and Red bands. It calculates $X_i = \ln(Value - L_{\infty})$.
Build Vector Y: Vector Y contains the known depths $z$ from the CSV.
Regression: A multiple linear regression is performed (function multipleLinearRegression) to find the coefficient vector $\beta$ (containing $a_0, a_1, a_2$).
Application: The found coefficients are applied to every pixel in the image to generate the bathymetric map.