A Likelihood-based Triangulation Method for Uncertainties in Through-Water Depth Mapping

1. Introduction

The coastal environment is of major interest to a broad spectrum of institutional structures as well as scientific communities and industrial entities. This is largely due to the economic and social stakes (ex. harbour, urban areas), strategic issues and the commercial activities (ex.aquaculture) concentrated in these areas [1,2]. Bathymetry is vital for coastal engineering studies as it provides information on the morphology of the seabed. Indeed, bathymetry-derived products are crucial for most maritime applications: navigation safety (navigation charts compilation), harbor development (pre/post-dredging surveys), and monitoring aquatic resources. On the technological side, various innovations emerged in the last decades leading to the development of sensors dedicated to bathymetric mapping [3]. Innovations in bathymetry include optical alternatives for shallower waters, as traditional acoustic solutions are less efficient and require greater investments.

Table 1. Parameter Symbols and Descriptions.

Symbol	Description
$\mathit{\theta }$	The parameter underlying the data (multidimensional)
$\tilde{\mathit{\theta }}$	True value of the parameter $\mathit{\theta }$
$\psi$	Parameter of interest (scalar)
$\mathit{\omega }$	Nuisance vector parameter (multidimensional)
$\mathbf{y}$	Observed data (multidimensional)
$\ell (\mathit{\theta };\mathbf{y})$	Likelihood function
$\ell (\mathit{\theta };\mathbf{y})$	Log-likelihood function
${\ell }_p(\mathit{\theta };\mathbf{y})$	Profile Log-likelihood function
${\ell }_{\theta \theta }\left[(\mathit{\theta };\mathbf{y})\right]$	Second derivative of the log likelihood with respect to $\mathit{\theta }$
$\widehat{\mathit{\theta }}\left(\mathbf{y}\right)$	Maximum likelihood estimate of $\mathit{\theta }$ given $\mathbf{y}$
$j_p(\psi ;\mathbf{y})$	Partial Observed Fisher information
$t_o(\psi ;\mathbf{y})$	Wald statistic (observed)
$r(\psi ;\mathbf{y})$	Signed root likelihood ratio statistic
$\alpha$, $1-\alpha$	Significance level, Confidence level (CL)
${\mathit{v}}_{\mathit{s}}$	Interior orientation vector
${\mathit{v}}_{\mathit{f}}$	Feature vector
${\mathit{x}}_{\mathit{f}}$	Feature point position
${\mathit{x}}_{\mathit{i}}$	Incidence point position
${\mathit{v}}_{\mathit{i}}$	Line-of-sight vector
${\mathit{q}}_{\mathit{i}}$	Line-of-sight quaternion
h	Water-air interface height parameter
$n_{\mathrm{water}},n_{\mathrm{air}}$	Refraction indices for water and air medium.
${\xi }_i$,${\xi }_r$	Incidence angle and refracted angle
${\xi }_b$	Backward refraction angle
${\mathit{v}}_{\mathit{b}}$	Backwardline-of-sight vector
$h-{\left[{\mathit{x}}_f\right]}_z$	Water Column Depth (WCD)
${\mathit{n}}_{\mathit{z}}$	Nadir vector
${\mathit{x}}_{\mathit{s}}$	Parameter of camera position vector
${\mathit{x}}_s$	Measured camera position vector
${\Sigma }_{\mathit{s}}$	Variance-covariance matrix of camera position
$\mathbf{q}$	Measured line-of-sight quaternion
$\mathit{M}$	Bingham orientation parameter
$\mathit{C}$	Bingham concentration parameter

Table 1. Parameter Symbols and Descriptions.

Symbol	Description
$\mathit{\theta }$	The parameter underlying the data (multidimensional)
$\tilde{\mathit{\theta }}$	True value of the parameter $\mathit{\theta }$
$\psi$	Parameter of interest (scalar)
$\mathit{\omega }$	Nuisance vector parameter (multidimensional)
$\mathbf{y}$	Observed data (multidimensional)
$\ell (\mathit{\theta };\mathbf{y})$	Likelihood function
$\ell (\mathit{\theta };\mathbf{y})$	Log-likelihood function
${\ell }_p(\mathit{\theta };\mathbf{y})$	Profile Log-likelihood function
${\ell }_{\theta \theta }\left[(\mathit{\theta };\mathbf{y})\right]$	Second derivative of the log likelihood with respect to $\mathit{\theta }$
$\widehat{\mathit{\theta }}\left(\mathbf{y}\right)$	Maximum likelihood estimate of $\mathit{\theta }$ given $\mathbf{y}$
$j_p(\psi ;\mathbf{y})$	Partial Observed Fisher information
$t_o(\psi ;\mathbf{y})$	Wald statistic (observed)
$r(\psi ;\mathbf{y})$	Signed root likelihood ratio statistic
$\alpha$, $1-\alpha$	Significance level, Confidence level (CL)
${\mathit{v}}_{\mathit{s}}$	Interior orientation vector
${\mathit{v}}_{\mathit{f}}$	Feature vector
${\mathit{x}}_{\mathit{f}}$	Feature point position
${\mathit{x}}_{\mathit{i}}$	Incidence point position
${\mathit{v}}_{\mathit{i}}$	Line-of-sight vector
${\mathit{q}}_{\mathit{i}}$	Line-of-sight quaternion
h	Water-air interface height parameter
$n_{\mathrm{water}},n_{\mathrm{air}}$	Refraction indices for water and air medium.
${\xi }_i$,${\xi }_r$	Incidence angle and refracted angle
${\xi }_b$	Backward refraction angle
${\mathit{v}}_{\mathit{b}}$	Backwardline-of-sight vector
$h-{\left[{\mathit{x}}_f\right]}_z$	Water Column Depth (WCD)
${\mathit{n}}_{\mathit{z}}$	Nadir vector
${\mathit{x}}_{\mathit{s}}$	Parameter of camera position vector
${\mathit{x}}_s$	Measured camera position vector
${\Sigma }_{\mathit{s}}$	Variance-covariance matrix of camera position
$\mathbf{q}$	Measured line-of-sight quaternion
$\mathit{M}$	Bingham orientation parameter
$\mathit{C}$	Bingham concentration parameter

Technological innovations over the last few decades have led to the development of various cameras dedicated to bathymetric mapping. A significant amount of research involving spectral imagery (hyperspectral and multispectral) has shown the ability to retrieve water biophysical properties such as: chlorophyll-A concentration, suspended particulate matter, benthic habitats in addition to the WCD. There are two major approaches for retrieving WCD from spectral imagery: radiometric and geometric. Radiometric approaches require spectral measurements as input for providing an estimation of WCD at a pixel level [4,5,6,7]. They have been mostly applied on satellite imagery where they contributed to the establishment of the Satellite Derived Bathymetry (SDB). On the other hand, geometric approaches allow to obtain 3D reconstruction from images by triangulating a set of pixels that are located on more than one image (homologous feature points). Recent studies demonstrated the potential of various photogrammetric (i.e. geometric) techniques for mapping shallow water depths from satellite imagery [8,9,10] and airborne imagery [11,12]. [10] reported similar accuracy to established SDB methods and argue that through water photogrammetry is a useful option in circumstances where radiometric methods are not applicable due to heterogeneous environments, inaccurate atmospheric correction or unavailable in-situ depth measurements.

As the navigation safety of vessels relies on WCD accuracy, it is of utmost importance to quantify the uncertainty associated with the derived product. Studies on bathymetry estimation using geometric approaches tend to focus on the estimation accuracy, rather than the assessment of uncertainty, despite the fact that there are many sources of uncertainty. Extending the principles involved in land photogrammetry to through-water photogrammetry presents the challenge imposed by integrating the refraction into the collinearity equations. Furthermore, it relies on the assumption of negligible local errors that is suitable for applying uncertainty propagation as recommended in [13]. As a result, there is a knowledge gap in the literature in terms of methods for estimating the uncertainty of the WCD calculated using geometric approaches in shallow coastal areas.

In this work, we introduce a novel methodology for assessing the WCD uncertainty computed through stereo-photogrammetric airborne imagery. Our approach differs significantly from collinearity-based triangulation methods, particularly in its suitability for coastal area surveys. Coastal environments present unique challenges, including limited access, dynamic landscapes, safety concerns, and environmental sensitivities, which complicate the establishment of Ground Control Points (GCP). These points are crucial for aligning imagery and refining the exterior orientation parameters. In coastal areas where GCP are scarce, relying on direct georeferencing with onboard Inertial Navigation Systems (INS) may offer limited accuracy compared to GCP-assisted land photogrammetry. This reduces the ability to correct exterior orientation parameters, potentially compromising the precision of the WCD estimate. Our methodology, therefore, focuses on determining the true object coordinates based on the observed camera pose, assuming that the interior orientation parameters are fixed and well-characterized. Additionally, the likelihood-based framework broadens the scope of traditional variance-covariance approaches by offering a more extensive uncertainty analysis, especially in complex scenarios such as small sample sizes where traditional approaches might underperform. Moreover, the use of controlled simulations instead of real datasets allows for a more objective assessment of our methodology’s robustness under various conditions, ensuring its validity and practical applicability in real world scenarios.

The structure of this paper is organized as follows: The first section provides an overview of existing approaches for uncertainty characterization in stereo-photogrammetry and highlights their limitations in airborne marine applications. The second section discusses the probabilistic framework employed for evaluating stereo-triangulated WCD uncertainties. The third section presents experimental results and discusses the performance of our proposed methodology under simulated scenarios. Finally, the conclusion summarizes the key findings, underscores the contributions and the potential of our likelihood-based triangulation approach, and discusses future research directions in this rapidly evolving field.

2. Related Works

This section synthesizes current methodologies for geometric WCD inference with a focus on uncertainty evaluation in through-water photogrammetry—an emergent field with a lack in comprehensive studies on uncertainty despite its importance alongside empirical validations.

Rational Polynomial Coefficients (RPC) are typically used in satellite pushbroom datasets for triangulating feature points and deriving bathymetric estimations in shallow waters [8,10]. Dolloff et al. have demonstrated the use of the variance-covariance framework to derive triangulation uncertainties using RPC data in land contexts [14,15]. However, in through-water contexts, the uncertainty inherent in RPC-based methods presents significant challenges. In principle, data on RPC uncertainty is essential to apply the variance-covariance framework for deriving triangulation uncertainties. In practice, most procedures for quantifying RPC uncertainties rely on GCP field measurements, which are challenging to obtain in coastal areas [16]. Furthermore, the need for refraction correction when using RPC-based triangulation introduces additional complexities. Although it is conceptually possible to implement the variance-covariance framework on top of RPC-based triangulation that includes refraction adjustments, this process adds a layer of complication to the uncertainty analysis in these settings. This lack of direct applicability complicates the uncertainty analysis for through-water photogrammetry, emphasizing the need for tailored methodologies that can accurately incorporate these refraction adjustments.

Structure from Motion (SfM) and Multi View Stereo (MVS) methods, while well-established in land photogrammetry, have been adapted for through-water use to account for the unique challenges posed by the aquatic environment [17,18,19]. In the terrestrial context, these methods essentially employ the variance-covariance propagation with the associated LS estimator to derive uncertainties [20]. However, when adapted to aquatic contexts, these methods deploy adaptations of the collinearity equations to accommodate the refraction effect. These required adaptations often involve adjustments either in the image space [12,21] or the object space [11,22], aiming to estimate refraction effects accurately. Such adjustments are crucial because adapting the variance-covariance approaches to these modified settings is not straightforward. Although no studies were found detailing a direct application of variance-covariance to these adaptations, the required iterative processes go beyond classical least squares methods, complicating the evaluation of the uncertainties. Furthermore, implicitly compensating for refraction in these approaches renders the physical interpretability of the results more challenging as the refraction effect is not explicitly accounted for. Besides the assumption of a flat water surface, these methods are conceived for scenarios where the camera is vertically oriented to the surface [23,24]. This suggests their sub-optimality when drone imagery, usually characterized by low stability, is being used to create a 3D model of the seabed.

[23] contrasts corrective methods, i.e. adaptations of RPC and SfM-MVS for through-water use, with ray tracing approaches, highlighting their limitations in handling complex optical paths effectively. This comparison emphasizes the need for tailored methodologies that can more effectively incorporate complex refraction adjustments while maintaining robust uncertainty analysis in through-water photogrammetry. Ray tracing methods, explicitly modeling refraction, offer a rigorous alternative for through-water triangulation and geometric uncertainties evaluation. Significantly, the explicit modeling employed by these methods enables a segmented analysis of the uncertainty budget, which in turn allows for the extraction of meaningful insights regarding the quality of measurements. Furthermore, advancements indicate that optimizing the cost function in the object space enhances geometric accuracy compared to traditional methods that primarily focus on the image space. [23,25]. Despite these advancements, it is important to note that coastal through-water photogrammetry studies have not yet fully incorporated these techniques, largely due to the predominance of corrective methods dictated by common software and dataset conventions in practice.

The discussed methods underscore the need for novel inferential approaches that integrate both triangulation and uncertainty evaluation, specifically tailored for the unique challenges of coastal environments. Our proposed likelihood-based triangulation method addresses these challenges by offering an inclusive and extensible framework for evaluating uncertainties in WCD inferences, leveraging advancements in both theoretical and practical aspects of photogrammetry.

3. Methodology

This section delineates our approach for triangulating the positions of feature points in through-water stereo-photogrammetry and quantifying the uncertainties associated with WCD measurements. Our methodology mainly adopt the MLE framework, which is theoretically equivalent to the commonly used least squares under specific assumptions and thus is inclusive of traditional methods. Unlike the decoupled procedures of least squares followed by variance-covariance framework for uncertainty evaluation, the MLE framework integrates estimation and uncertainty quantification into a cohesive process. Its statistical properties, robust in large samples supported by first-order theory, make it a comprehensive joint framework for both estimation and uncertainty analysis. Following this, the first section 3.1 addresses the proposed Likelihood triangulation and the associated statistical modeling assumptions. This allows for establishing the likelihood function for through-water triangulation as well as the different geometric modeling aspects of the proposed approach. The second section addresses the statistical evaluation of the uncertainties, employing first-order statistical tests for parameter uncertainty.

3.1. Proposed Likelihood Triangulation

Classical triangulation in photogrammetry typically utilizes collinearity equations to relate the object coordinates ${\mathit{x}}_{\mathit{f}}$ of a feature point to its image coordinates ${\mathit{x}}_{\mathit{p}}$, and the camera’s interior and exterior orientations, whose parameters can be denoted as ${\theta }_{\mathit{io}}$ and ${\theta }_{\mathit{eo}}$ respectively.

As depicted in Figure 1, the primary goal in this context is to minimize reprojection error between the observed and modeled feature image coordinates, denoted as ${\mathit{x}}_{\mathit{p}}$. This is formally achieved through the following LS minimization:

$$min{\left({\mathit{x}}_p-{\mathit{x}}_{\mathit{p}}({\mathit{x}}_{\mathit{f}},{\theta }_{\mathit{io}},{\theta }_{\mathit{eo}})\right)}^2$$

(1)

Here, ${\mathit{x}}_p$ denotes the observed image coordinates of the feature point. As a consequence of the equivalence between the least squares and the MLE when assuming additive Gaussian errors, this optimization effectively seeks to maximize the likelihood:

$$L\left({\mathit{x}}_{\mathit{p}}({\mathit{x}}_{\mathit{f}},{\theta }_{\mathit{io}},{\theta }_{\mathit{eo}});{\mathbf{x}}_p\right)$$

(2)

This reprojection likelihood function describes the probability of observing the feature point image coordinates given the parameters ${\mathit{x}}_{\mathit{f}}$, ${\theta }_{\mathit{eo}}$, ${\theta }_{\mathit{io}}$. To achieve accurate triangulation, calibration is critical to optimize interior orientation parameters such as the focal length across surveys, while exterior orientations are often refined through spatial resection using GCP to enhance accuracy. Triangulation methods often rely on optimal values of ${\theta }_{\mathit{io}}$ and ${\theta }_{\mathit{eo}}$, disregarding the integration of the uncertainty of interior and exterior orientations into the estimation of the triangulated position of the feature points. This limitation is particularly critical in coastal surveys, where access to GCP is often restricted. It becomes even more crucial for pushbroom cameras, which must resolve or refine exterior orientation per swath rather than per image. Additionally, the minimal texture and feature distinctiveness of seabeds in coastal waters may complicate pairing of features across overlapping images, further reducing confidence in the exterior orientation.

Figure 1. Illustration of the process of minimizing reprojection error in stereo-photogrammetry. Diagram of camera setup showing the Field Of View (FoV), focal length, and the Instantaneous Field Of View (IFoV). The internal vector ${\mathit{v}}_{\mathit{s}}$ and feature image coordinates ${\mathit{x}}_{\mathit{p}}$ are highlighted, with a color bar gradient indicating the reprojection likelihood magnitude.

Figure 1. Illustration of the process of minimizing reprojection error in stereo-photogrammetry. Diagram of camera setup showing the Field Of View (FoV), focal length, and the Instantaneous Field Of View (IFoV). The internal vector ${\mathit{v}}_{\mathit{s}}$ and feature image coordinates ${\mathit{x}}_{\mathit{p}}$ are highlighted, with a color bar gradient indicating the reprojection likelihood magnitude.

To address these challenges, our methodology emphasizes probabilistic modeling of camera pose, focusing particularly on its impact on the spatial uncertainty of the line-of-sight from the camera to the incidence point ${\mathit{v}}_{\mathit{i}}$. Figure 2 illustrates this likelihood function which can be employed in order to estimate the feature point position from camera pose data:

$$L({\mathit{v}}_{\mathit{i}}(h,{\mathit{x}}_{\mathit{f}},{\mathit{x}}_{\mathit{s}});\mathbf{y}|{\mathit{v}}_{\mathit{s}})$$

(3)

Here, $\mathbf{y}$ represents the measured camera pose data, ${\mathit{v}}_{\mathit{s}}$ is a unit vector representing the feature image coordinates in the sensor frame, h is the Water Air Interface (WAI) height and ${\mathit{x}}_{\mathit{s}}$ is the modeled camera source position. Since the influence of interior orientation is disregarded in the proposed approach, the internal vector ${\mathit{v}}_{\mathit{s}}$ which describes the feature coordinates in the sensor frame is considered fixed and known posterior to feature detection and pairing. This allows to isolate the influence of the camera pose noise on the line-of-sight ${\mathit{v}}_{\mathit{i}}$ and establish its spatial likelihood as illustrated in Figure 2.

Figure 2. Line-of-sight representation in through-water photogrammetry. Representation of the camera’s line-of-sight ${\mathit{v}}_{\mathit{i}}$, demonstrating its path through the water column to the seabed feature position ${\mathit{x}}_{\mathit{f}}$. Color bar gradient illustrates the magnitude of the line-of-sight likelihood.

Figure 2. Line-of-sight representation in through-water photogrammetry. Representation of the camera’s line-of-sight ${\mathit{v}}_{\mathit{i}}$, demonstrating its path through the water column to the seabed feature position ${\mathit{x}}_{\mathit{f}}$. Color bar gradient illustrates the magnitude of the line-of-sight likelihood.

In this study, we focus on bathymetric measurements, particularly examining both the WCD and WAI height presented in Figure 3. The WCD describes the geometry of underwater environments, while the WAI marks the boundary between water and air, and its ellipsoidal height h is used for surface refraction modeling.

3.1.1. Line-of-Sight Modeling

For each environmental geometric configuration consisting of the feature position ${\mathit{x}}_{\mathit{f}}$ and WAI height h, the objective is to model the line-of-sight vector ${\mathit{v}}_{\mathit{i}}$ using the parameters of the feature object coordinate ${\mathit{x}}_{\mathit{f}}$, the camera source position ${\mathit{x}}_{\mathit{s}}$, the horizontal water-air interface height h, and the nadir ${\mathit{n}}_{\mathit{z}}$. This allows to account for refraction effect at the WAI through a strict geometric representation.

The steps involved in modeling the line of sight in the presence of a horizontal WAI are as follows:

1.

The incidence angle ${\xi }_i$,the angle between the feature vector ${\mathit{v}}_{\mathit{f}}={\displaystyle \frac{{\mathit{x}}_{\mathit{f}}-{\mathit{x}}_{\mathit{s}}}{|{\mathit{x}}_{\mathit{f}}-{\mathit{x}}_{\mathit{s}}|}}$ and the nadir ${\mathit{n}}_z$, is calculated as:

$${\xi }_i={cos}^{-1}({\mathit{v}}_{\mathit{f}}·{\mathit{n}}_{\mathit{z}})$$

(4)

2.

Using Snell’s Law, the refraction angle ${\xi }_r$ when transitioning from water to air is:

$${\xi }_r={sin}^{-1}(nsin\left({\xi }_i\right))$$

(5)

where n is the refractive index ratio, fixed at 1.33, representing the ratio of the speed of light in air to that in water. The incidence and refraction angles calculated in steps 1 and 2 are computed as if there was no refraction. These angles do not represent the actual incidence and refraction angles in the presence of refraction, but rather serve as proxies to determine the backward rotation necessary to model the refraction effect in the subsequent steps.

3.

The adjusted vector ${\mathit{v}}_{\mathit{b}}$ to refraction, a function of ${\mathit{v}}_{\mathit{f}}$, ${\mathit{n}}_{\mathit{z}}$, and the refractive index n, is computed by applying a rotation ${\mathit{R}}_b\left({\xi }_b\right)$ in the plane ${\mathit{n}}_z\times {\mathit{v}}_{\mathit{i}}$, witg angle ${\xi }_b={\xi }_i-{\xi }_r$:

$${\mathit{v}}_{\mathit{b}}={\mathit{R}}_b\left({\xi }_b\right)·(-{\mathit{v}}_{\mathit{f}})$$

(6)

4.

The incidence point ${\mathit{x}}_{\mathit{i}}$ on the interface, where ${\mathit{v}}_{\mathit{b}}$ intercepts the water-air interface, is calculated as:

$${\mathit{x}}_{\mathit{i}}={\mathit{x}}_{\mathit{f}}+\left({\displaystyle \frac{h-{\left[{\mathit{x}}_f\right]}_z}{{\left[{\mathit{v}}_{\mathit{b}}\right]}_z}}\right){\mathit{v}}_{\mathit{b}}$$

(7)

where ${\left[{\mathit{x}}_f\right]}_z$ and ${\left[{\mathit{v}}_b\right]}_z$ are the vertical componenets of ${\mathit{x}}_{\mathit{f}}$ and ${\mathit{v}}_{\mathit{b}}$ respectively.

5.

Finally, the line-of-sight vector ${\mathit{v}}_{\mathit{i}}$, the normalized vector from ${\mathit{x}}_{\mathit{s}}$ to ${\mathit{x}}_{\mathit{i}}$, is defined as:

$${\mathit{v}}_{\mathit{i}}={\displaystyle \frac{{\mathit{x}}_{\mathit{i}}-{\mathit{x}}_{\mathit{s}}}{|{\mathit{x}}_{\mathit{i}}-{\mathit{x}}_{\mathit{s}}|}}$$

(8)

Figure 3 illustrates the flow of Equations (4)–(8) which collectively define the model for the line-of-sight vector ${\mathit{v}}_{\mathit{i}}$. This provides a direct relationship between the optical path and the geometric parameters h, ${\mathit{x}}_{\mathit{f}}$, and ${\mathit{x}}_{\mathit{s}}$ and ensures a strict modeling of through-water refractive geometry. This model will now be used to address the influence of the camera’s position and orientation on the spatial uncertainty of the line-of-sight vector.

3.1.2. Camera Pose Statistical Model

Given the geometric model for the line-of-sight vector ${\mathit{v}}_{\mathit{i}}$, we can directly connect these calculations to the triangulation process, pivotal for accurate underwater reconstructions. Considering a known interior orientation, the line-of-sight vector ${\mathit{v}}_{\mathit{i}}$ depends only on the camera pose.

3.1.2.a Camera Pose Data

The position of the camera’s optical center is denoted ${\mathit{x}}_{\mathit{s}}$ and its orientation is modeled using the unit quaternion ${\mathit{q}}_{\mathit{i}}$ representing the minimum rotation between the internal vector ${\mathit{v}}_{\mathit{s}}$ and the line-of-sight ${\mathit{v}}_{\mathit{i}}$ pointing towards the incidence point as illustrated in Figure 4. Unit quaternions are selected for representing the camera’s orientation to avoid the limitations associated with Euler angles, such as gimbal lock and the complexity of interpolation. Unit quaternions also provide a compact, non-singular representation for rotations and are particularly effective for continuous orientation tracking in three-dimensional space.

The correspondence between the minimum rotation unit quaternions and the lines-of-sights is unique up to sign meaning that both ${\mathit{q}}_{\mathit{i}}$ and $-{\mathit{q}}_{\mathit{i}}$ represent the same orientation. The unit quaternion ${\mathit{q}}_{\mathit{i}}$ defining the minimal rotation from ${\mathit{v}}_{\mathit{s}}$ to ${\mathit{v}}_{\mathit{i}}$ can be defined as:

$$\begin{array}{c}{q}_w=|{\mathit{v}}_{\mathit{s}}\left|\right|{\mathit{v}}_{\mathit{i}}|+{\mathit{v}}_{\mathit{s}}·{\mathit{v}}_{\mathit{i}}\end{array}$$

(9)

$$\begin{array}{c}\left[\begin{array}{ccc}{q}_x& q_y& q_z\end{array}\right]={\mathit{v}}_{\mathit{s}}\times {\mathit{v}}_{\mathit{i}}\end{array}$$

(10)

$$\begin{array}{c}{\mathit{q}}_{\mathit{i}}={\displaystyle \frac{1}{\sqrt{q_x^2+q_y^2+q_z^2+q_w^2}}}\left[\begin{array}{cccc}{q}_x& q_y& q_z& q_w\end{array}\right]\end{array}$$

(11)

$|{\mathit{v}}_{\mathit{s}}|$ and $|{\mathit{v}}_{\mathit{i}}|$ are the magnitudes of the internal vector and the line-of-sight , respectively and ${\mathit{v}}_{\mathit{s}}·{\mathit{v}}_{\mathit{i}}$ is the dot product of both vectors. ${\mathit{v}}_{\mathit{s}}\times {\mathit{v}}_{\mathit{i}}$ is the cross product yielding a vector that is perpendicular to both input vectors, which represents quaternion components $q_x$, $q_y$, and $q_z$ as the axis rotation component. ${\mathit{q}}_{\mathit{i}}$ is normalized to have unit magnitude.

3.1.2.b Statistical Model

The position of the camera ${\mathit{x}}_s$ is modeled using a Gaussian distribution due to its properties in reflecting Euclidean space errors:

$$f_{X_s}({\mathit{x}}_s;{\mathit{x}}_{\mathit{s}},{\Sigma }_s)={\displaystyle \frac{1}{\sqrt{{\left(2\pi \right)}^{3}det\left({\Sigma }_s\right)}}}exp\left(-{\displaystyle \frac{1}{2}}{({\mathit{x}}_s-{\mathit{x}}_{\mathit{s}})}^T{\Sigma }_s^{-1}({\mathit{x}}_s-{\mathit{x}}_{\mathit{s}})\right)$$

(12)

Here, ${\mathit{x}}_{\mathit{s}}$ is the measured position vector of the sensor and ${\mathit{x}}_{\mathit{s}}$ and ${\Sigma }_{\mathit{s}}$ are parameters of the sensor position and its variance-covariance matrix, respectively.

The camera orientation determined by a unit quaternion $\mathbf{q}$ is modeled using the Bingham distribution, which is specifically designed for data constrained to the surface of a hypersphere, like unit quaternions. This distribution is parameterized by an orientation matrix $\mathit{M}$ and a concentration matrix $\mathit{C}$, which adjust the distribution’s shape based on the degree of certainty in the camera’s orientation:

$$f_Q(\mathbf{q};\mathit{M},\mathit{C})={\displaystyle \frac{1}{F\left(\mathit{C}\right)}}exp\left({\mathbf{q}}^T\mathit{M}\mathit{C}{\mathit{M}}^T\mathbf{q}\right)$$

(13)

Here, $\mathbf{q}$ is the random attitude measurement as the quaternion that represents the observed orientation of the sensor. In the Bingham distribution, $\mathit{M}\in {\mathbb{R}}^{4\times 4}$ is an orthogonal matrix and $\mathit{C}$ is a diagonal matrix of size $4\times 4$ representing the orientation and concentration parameters, respectively.

The concentration parameter $\mathit{C}$ determines the spread or the shape of the distribution, whereas the orientation parameter $\mathit{M}$ represents the location of the distribution, also called the orientation. The diagonal elements of the concentration matrix $(C_1\le C_2\le C_3\le C_4=0)$ control the degree of anisotropy of the distribution. Due to the constraint of unit norm on the quaternions, one of the diagonal concentration values must be zero. The fourth diagonal element $C_4$ is chosen while other conventions may use the first element. An isotropic distribution has equal concentration along all three axes of rotation, resulting in a uniform distribution of orientations $C_1=C_2=C_3$ (see Figure 5). In an anisotropic distribution, the diagonal elements of $\mathit{C}$ are not equal, and the distribution is stretched along one or more axes of rotation. $F\left(\mathit{C}\right)$ is a normalization constant that depends on the concentration of the density.

At this stage, it is critical to link the line-of-sight ${\mathit{v}}_i$ to the Bingham distribution parameters to accurately formulate its likelihood. The orientation matrix $\mathit{M}$, where its last column is defined by the quaternion ${\mathit{q}}_i({\mathit{v}}_i,{\mathit{v}}_{\mathit{s}})$, plays a pivotal role. This quaternion aligns the sensor’s internal frame, described by vector ${\mathit{v}}_{\mathit{s}}$, with the line-of-sight vector ${\mathit{v}}_i$. The remaining columns of $\mathit{M}$ in isotropic conditions can be filled with an arbitrary orthogonal complement using methods such as the Gram-Schmidt process or singular value decomposition of ${\mathit{q}}_i{\mathit{q}}_i^T$. Orientation data from (INS), characterized by pronounced anisotropic uncertainties, with yaw being more variable than roll or pitch, necessitates sophisticated configuration of the Bingham distribution’s parameters, $\mathit{M}$ and $\mathit{C}$. Such data typically arrives in Eulerian format, with mean and covariance values computed through spatial resection techniques or from outputs of a Kalman filter, and is foundational in determining these parameters.

One can consider dynamically recomputing $\mathit{M}$ and $\mathit{C}$ are for each line-of-sight by directly sampling from the Eulerian noise associated with that specific orientation on the sphere. Although this appraoch is capable of capturing true variability in sensor orientation ensuring fidelity, it also requires sampling quaternion data points from the Eulerian data distribution. However, its computational intensity limits scalability across large datasets, particularly when considering optimiation tasks. For practicality in extensive datasets, we use this simpler approach. ${\mathit{M}}_a$ and ${\mathit{C}}_a$ are predetermined from a representative average of the orientation data by resampling from the Eulerian attitude noise. Additionnally, this static setup efficiently approximates the dynamic method under conditions of minimal attitude noise. The orientation density approximation can be formulated as:

$$f_Q(\mathbf{q};{\mathit{v}}_{\mathit{i}},{\mathit{v}}_{\mathit{s}})={\displaystyle \frac{1}{\mathit{F}\left({\mathit{C}}_{\mathit{a}}\right)}}exp\left({\mathit{q}}_i^T({\mathit{v}}_i,{\mathit{v}}_{\mathit{s}}){\mathit{M}}_a{\mathit{C}}_a{\mathit{M}}_a^T{\mathit{q}}_i({\mathit{v}}_i,{\mathit{v}}_{\mathit{s}})\right),$$

(14)

The matrices ${\mathit{M}}_a$, ${\mathit{C}}_a$, and the normalization factor $F\left({\mathit{C}}_a\right)$ are derived from the camera orientation data realization $\mathit{q}$ and a fixed Euler attitude covariance using the maximum likelihood estimation method described in [26] with 10000 randomly sampled quaternions. This ensures that the fixed parameters effectively model the orientation uncertainties, streamlining computation while maintaining analytical rigor in the evaluation of WCD uncertainties in through-water triangulation.

Having rigorously defined the statistical models for the camera’s position and orientation using Gaussian and Bingham distributions, we have established a robust framework to represent the stochastic nature of the camera’s pose. The line-of-sight vector ${\mathit{v}}_{\mathit{i}}$, as modeled in Section 3.1.1, depends on the horizontal water-air interface height h, the camera source position ${\mathit{x}}_{\mathit{s}}$, and the feature point coordinates ${\mathit{x}}_{\mathit{f}}$. Hence the position and the orientation distributions combine into a joint probability density function for the camera pose parameterized as follows:

$$f_Y(\mathbf{y};h,{\mathit{x}}_{\mathit{f}},{\mathit{x}}_{\mathit{s}})=f_{X_s}({\mathit{x}}_{\mathit{s}};{\mathit{x}}_{\mathit{s}},{\Sigma }_s)\times f_Q(\mathbf{q};{\mathit{v}}_{\mathit{i}},{\mathit{v}}_{\mathit{s}}),$$

(15)

where $f_Y$ encapsulates the statistical model of the camera’s position and orientation $({\mathbf{x}}_s,\mathbf{q})$ as a function of the line-of-sight parameters $h,{\mathit{x}}_{\mathit{f}},{\mathit{x}}_{\mathit{s}}$ since the internal vector ${\mathit{v}}_{\mathit{s}}$ is fixed. This joint density function provides a quantitative measure of uncertainty and forms the link to the line-of-sight vectors employed for through-water triangulation.

3.1.3. MLE Based Triangulation

Leveraging the presented camera pose statistical model, we construct the likelihood function that accounts for both positional and rotational uncertainties in the camera setup. The likelihood function L for a single feature point observed from a single camera can be expressed as:

$$L(h,{\mathit{x}}_{\mathit{f}},{\mathit{x}}_{\mathit{s}};\mathbf{y})=f_Y(\mathbf{y};h,{\mathit{x}}_{\mathit{f}},{\mathit{x}}_{\mathit{s}})$$

(16)

Assuming independent camera poses, the likelihood function L for a single feature point observed from a multiple cameras can be expressed as:

$$L(h,{\mathit{x}}_{\mathit{f}},{\mathit{x}}_{\mathit{s}};\mathbf{y})=\prod _{m=1}^M{f}_Y({\mathbf{y}}^m;h,{\mathit{x}}_{\mathit{f}}^m,{\mathit{x}}_{\mathit{s}}^m)$$

(17)

Disregarding the constant terms which do not involve the parameters, mainly the normalisation constants $\sqrt{{\left(2\pi \right)}^{3}det\left({\Sigma }_s\right)}$ and $\mathit{F}\left({\mathit{C}}_{\mathit{a}}\right)$, we can write:

$$L(h,{\mathit{x}}_{\mathit{f}},{\mathit{x}}_{\mathit{s}};\mathbf{y})\propto exp\left(-\sum _{m=1}^M\left[{({\mathbf{x}}_s^m-{\mathit{x}}_{\mathit{s}}^m)}^T{\Sigma }_{s,m}^{-1}({\mathbf{x}}_s^m-{\mathit{x}}_{\mathit{s}}^m)+{\mathit{q}}_i^T{\mathit{M}}_{a,m}{\mathit{C}}_{a,m}{\mathit{M}}_{a,m}^T{\mathit{q}}_{\mathit{i}}\right]\right)$$

(18)

The logarithm of this likelihood function, denoted as ℓ, can be written, ignoring the constant terms, as a single index m summation with two quadratic forms:

$$\ell (h,{\mathit{x}}_{\mathit{f}},{\mathit{x}}_{\mathit{s}}:\mathbf{y})=-\sum _{m=1}^M\left[{({\mathbf{x}}_s^m-{\mathit{x}}_{\mathit{s}}^m)}^T{\Sigma }_{s,m}^{-1}({\mathbf{x}}_s^m-{\mathit{x}}_{\mathit{s}}^m)+{\mathit{q}}_{i,m}^T{\mathit{M}}_{a,m}{\mathit{C}}_{a,m}{\mathit{M}}_{a,m}^T{\mathit{q}}_{i,m}\right]$$

(19)

For multiple feature points observed from multiple sensors, the logarithm of this likelihood function can be expressed as a double summation over k and m, where K represents the number of feature points and M represents the number of sensors:

$$\ell (h,{\mathit{x}}_{\mathit{f}},{\mathit{x}}_{\mathit{s}};\mathbf{y})=-\sum _{m=1}^M\sum _{k=1}^K\left[{({\mathbf{x}}_s^m-{\mathit{x}}_{\mathit{s}}^m)}^T{\Sigma }_{s,m}^{-1}({\mathbf{x}}_s^m-{\mathit{x}}_{\mathit{s}}^m)+{\mathit{q}}_{i,k,m}^T{\mathit{M}}_{a,m}{\mathit{C}}_{a,m}{\mathit{M}}_{a,m}^T{\mathit{q}}_{i,k,m}\right]$$

(20)

In a typical case, with two camera views observing a single feature point, the number of parameters to estimate is 9 including 3 for ${\mathit{x}}_{\mathit{f}}$, 1 for h, and 6 for both camera positions ${\mathit{x}}_{\mathit{s}}=({\mathit{x}}_{\mathit{s}}^1,{\mathit{x}}_{\mathit{s}}^2)$. The MLE is then obtained by maximizing this function with respect to the parameters h, ${\mathit{x}}_{\mathit{f}}$, and ${\mathit{x}}_{\mathit{s}}$. For problems where we estimate the WAI height h, we utilize multi-feature multi-view likelihood to leverage more statistical information and investigate accuracy improvements.

3.2. Uncertainties Evaluation

Understanding the uncertainties associated with our estimations is crucial for robust decision-making and inference. We begin by exploring the concept of profile likelihood, which allows us to focus on parameters of interest while accounting for nuisance parameters. This leads to first-order statistical tests, where we introduce the signed likelihood ratio statistic and the Wald test statistic hypothesis testing and confidence intervals (CI) computation. Subsequently, we describe the procedure of assessing the performance of CI in terms of coverage probability, shedding light on the reliability of the first order statistical tests. Our investigation aims to provide insights into the uncertainty quantification process, essential for interpreting and utilizing the results of our triangulation approach effectively.

3.2.1. Profile Likelihood

The likelihood function for through-water triangulation depends on the parameters $h,{\mathit{x}}_{\mathit{f}},{\mathit{x}}_{\mathit{s}}$ which can be decomposed into two components $\mathit{\theta }=(\psi ,\mathit{\omega })$ depending on the investigated parameter. The concept of profile likelihood plays a critical role in statistical inference, particularly when focusing on a parameter of interest $\psi$ while accounting for nuisance parameters:

$$l_p(\psi ;\mathbf{y})=\underset{\omega }{max}l(\psi ,\mathit{\omega };\mathbf{y})$$

(21)

In essence, the profile likelihood ${\ell }_p(\psi ;\mathbf{y})$ simplifies the likelihood function by focusing on $\psi$ and adjusting the nuisance parameters $\mathit{\omega }$ accordingly. This is achieved by maximizing the likelihood with respect to the nuisance parameter, resulting in the constrained maximum likelihood ${\widehat{\mathit{\omega }}}_{\psi }$ for each value of $\psi$ (see Figure 6).

Notably, maximizing the profile likelihood is in accordance with the overall MLE of $\mathit{\theta }=(\psi ,\mathit{\omega })$:

$$\widehat{\mathit{\theta }}=arg\underset{\psi ,\omega }{max}\ell (\psi ,\mathit{\omega };\mathbf{y})=arg\underset{\psi }{max}{\ell }_p(\psi ;\mathbf{y})$$

(22)

3.2.2. First Order Statistical Tests

For scalar parameter $\psi$, first-order inference uses test statistics like the signed likelihood ratio statistic r and the Wald test statistic $t_o$ [27]:

$$\begin{array}{cc}\hfill r(\psi ;\mathbf{y})& =\mathrm{sign}(\widehat{\psi }-\psi )\sqrt{2({\ell }_p(\widehat{\psi };\mathbf{y})-\ell (\psi ;\mathbf{y}))},\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill t_o(\psi ;\mathbf{y})& =j_p{(\widehat{\psi };\mathbf{y})}^{-1/2}(\widehat{\psi }-\psi ),\hfill \end{array}$$

(24)

where $j_p(\widehat{\psi };\mathbf{y})=-{\displaystyle \frac{{\partial }^2{\ell }_p(\widehat{\psi };\mathbf{y})}{\partial {\psi }^2}}$ is the partial observed information. Under certain regularity conditions, these test statistics asymptotically follow a standard normal distribution under the null hypothesis:

$$st(\tilde{\psi };\mathbf{y})\stackrel{O\left(n^{-1/2}\right)}{\to }\mathcal{N}(0,1)\mathrm{for}st\in \{r,t_o\}.$$

(25)

It should be emphasized that scalar $\psi$ test statistics are in harmony with multidimensional inference. In those broader contexts, the test statistics converge to a chi-squared distribution, where the degrees of freedom match the dimensionality of the interest parameter. This ensures that focusing on a scalar parameter $\psi$ does not deviate from the traditional multidimensional inference context. Furthermore, the first-order arsenal includes the score and the Wald test based on the Expected Fisher Information the in addition to the test statistics investigated in this study. In particular, it is noteworthy to emphasize the equivalence between the Wald test based on the Expected Fisher Information and variance-covariance/LS under additive Gaussian errors (see Section A). This underscores the profile likelihood framework as an encompassing approach that integrates traditional variance-covariance practices. While this study excludes the use of this test statistic, as an analytical expression of the Expected Fisher Information matrix can not be readily obtained due to the use of the Bingham distribution, it is essential to highlight that the Expected and observed information Wald tests often behave similarly in large samples.

3.2.3. Evaluation of Confidence Interval Performance

The performance of CI for the interest parameter is evaluated in coverage probability. Coverage probability is defined as the frequency with which the true value of the parameter of interest, denoted $\tilde{\psi }$, falls within the estimated CI $\mathcal{A}\left(\mathbf{y}\right)$. Mathematically, coverage probability is expressed as:

$$P(\tilde{\mathit{\psi }}\in \mathcal{A}\left(\mathbf{y}\right))=1-\alpha$$

(26)

Here, $\alpha$ represents the significance level, set at 0.05 for a 95% confidence level (CL) vin this study. To construct CI, the following first-order test statistics are employed:

$$\begin{array}{cc}\hfill \mathrm{For}r:& \psi :\left|r(\psi ;\mathbf{y})\right|\le z_{\alpha /2},\hfill \end{array}$$

(27)

$$\begin{array}{cc}\hfill \mathrm{For}{t}_o:& \widehat{\psi }\pm z_{\alpha /2}\sqrt{j_p(\widehat{\psi };\mathbf{y})}.\hfill \end{array}$$

(28)

$z_{\alpha /2}$ is the critical value corresponding to the $\alpha /2$ tail of a standard normal distribution, generally approximated as 1.96 for the 95% CL. Figure 7 demonstrates how hypothesis testing and CI estimation are performed based on the profile likelihood for the 95% CL. In order to empirically validate the coverage probability of CI, we use Monte Carlo simulations to generate sample data $\mathbf{y}$ from the actual parameter values $\tilde{\mathit{\theta }}$. Specifically, we calculate the observed test statistic for each camera pose sample $\mathit{y}$ at the true $\tilde{\psi }$ value under the null hypothesis. An observed test statistic’s absolute value greater than $z_{\alpha /2}=1.96$ would indicate that the null hypothesis can be rejected at the $95\%$ CL. The ideal performance under simulation would be rejecting the true hypothesis 5% of the time, indicating that the confidence intervals derived from the investigated test statisic align with their expected theoretical coverage.

Unlike variance-covariance approaches, profiling can reveal the asymmetric behavior of CI in nonlinear setups, crucial for accurate uncertainty quantification. As highlighted in the Figure 7, Wald tests $t_o$ result by construction in a linear slope yielding symmetric CI around the MLE $\widehat{\psi }$. On the other hand, r effectively captures the asymetric behavior which can arise in non linear setups.

For the WCD uncertainty evaluation, we consider the constrained triangulation problem where h is fixed and known whereas the primary parameter of interest is the WCD, i.e. $h-{\left[x_f\right]}_z$. Additionally, we consider the case where the interest parameter is a systematic the WAI height h which is to be inferred from multiple features observations.

4. Results

In the this section, we present an overview of the main results obtained from our probabilistic modeling and likelihood-based inference approach. We first begin by describing the simulated experiments consisting in realistic pushbroom geometric configurations in coastal areas, used to validate the proposed methodology. The results are partitioned into two main sections: WCD inference and WAI height inference focusing on the uncertainties of these parameters.

4.1. Simulated Experiments

The experimental setup for this study is designed to evaluate the performance of the proposed method for pushbroom camera acquisition under various flight scenarios, encompassing both airborne and drone operations. The airborne scenario simulates an altitude of 2000 meters, reflecting typical conditions for aircraft-based remote sensing missions, while the drone scenario simulates a lower altitude of 120 meters, indicative of low-altitude drone surveys.

Figure 8 illustrates the camera positions (black squares) and their respective scanned FoV (blue and orange triangles) for the drone scenario. We use a FoV of 48°, a relatively extreme value to challenge our approach. To ensure comparability and consistency, the same base-height ratio determined for along-flight directions (determined by the 60% overlap constraint) is transposed to the single across-flight trajectory used for the $P_c$ scenario. This setup maintains a consistent geometric configuration across different orientations.

In the experimental setup, various viewing geometries were simulated to evaluate the performance of the proposed method. Table 2 presents these viewing geometries which are determined by feature points $P_l$, $P_r$ and $P_c$. Table 3 presents the simulated INS classes defining sensor pose quality metrics in terms of standard deviation.

For both the drone and the airborne experiments, three levels of camera pose quality were considered: "Fair", "Good", and "Excellent". The positional accuracy ranges from .5 meters for the Fair class to 0.05 meters for the "Good" and "Excellent" classes. In our simulations, the precision of the camera’s attitude is modeled with distinct values for roll/pitch and heading. Specifically, for the ’Fair’ class, we use 0.1 degrees for roll/pitch precision and 1 degree for heading precision. For the ’Good’ class, the precision is set at 0.01 degrees for roll/pitch and 0.1 degrees for heading. Lastly, the ’Excellent’ class is characterized by a uniform precision of 0.01 degrees for both roll/pitch and heading.

To assess the robustness and reliability of the proposed methodology, we perform Monte Carlo simulations by generating synthetic datasets around a true parameter value. A true attitude value of $\tilde{\mathit{a}}=(0,0,0)$ and a fixed covariance matrix ${\Sigma }_{\mathit{a}}$ were considered for the attitude while the true camera positions were simply set to the simulated locations. These simulations involve generating a large number of synthetic data sets represented by random orientation and concentration matrices ${\mathit{M}}_a$,${\mathit{C}}_a$, as well as random camera positions ${\mathbf{x}}_s$, each with different levels of noise. The results of the simulations are used to evaluate the performance of the methodology and to analyze the uncertainties of WCD and the WAI height in terms of 95% CI under different conditions. Concerning the numerical computation of the MLE, the Adam optimizer is used for the WAI height inference since the dimension of the parameter space is relatively large in this problem. The Trust Region optimizer available in the SciPy library is employed for the inference of WCD, chosen for its efficiency in handling problems of relatively small scales.

4.2. Water Column Depth Inference

In this section, we consider the WCD inference in the through-water context, assuming that the water air interface height is known and fixed at $h=0$.

4.2.1. WCD Uncertainties

Figure 9 provides a graphical representation of the normalized profile likelihood surfaces for the three points: $P_l$, $P_r$, and $P_c$, specifically for the drone case. Each plot corresponds to a planar cut along the yz-plane of the 3D multi-view likelihood profiled along the interest parameter ${\mathit{x}}_{\mathit{f}}$, displaying the confidence regions obtained using the likelihood ratio with a degree of freedom of 3. We observe that the uncertainties for $P_l$ are generally higher compared to $P_r$ and $P_c$ with a factor of 2 between $P_l$ and $P_r$ uncertainties. In the one-medium photogrammetry, as it is known that the vertical uncertainty is proportional to the base height ratio. Additionally, the point $P_c$, which is viewed by perpendicular lines, exhibits non vertical uncertainty ellipses mainly because the line-of-sights are not co-planar as suggested by the difference in camera contribution to the 3D likelihood in the yz-plane as highlighted by Figure 9 for point $P_c$. Analyzing the uncertainty ellipses within the experiment classes present in Table 3, our results indicate that the WCD uncertainty shape is not heavily influenced by the depth factor.

Focusing on the WCD as a parameter of interest, we can further examine its uncertainties based on $t_o$ statistic using Monte Carlo simulations (as shown in Figure 10). The plot in Figure 10 provides a comparative analysis of the WCD uncertainty in terms of 95% confidence interval widths along with error bars representing the standard deviation of these intervals.

Overall, the impact of depth on the uncertainty of WCD is less significant than other viewing geometries. However, there is a slight increase in uncertainty as the depth increases especially in the drone scenario , with 15% uncertainty increase from 1 meter to 20 meters depth vs 3% uncertainty increase in the airborne scenarios. This is mainly due to the fact that base height ratio decreases as the depth of the simulated feature point increases, which generates an elongation of uncertainty ellipses. In terms of camera pose quality, both the "Excellent" INS and "Good" INS scenarios demonstrate relatively low and similar uncertainty values (50 cm to 1 m in the drone case and 3 m to 6m in the airborne case), with similar trends observed across different points. This indicates that increasing yaw precision does not necessarily lead to a more precise WCD estimation given the defined viewing geometries. Yaw rotation has a minimal effect on vectors which are close to the nadir suggesting that the internal vectors ${\mathit{v}}_{\mathit{s}}$ are not influenced by yaw precision even with a high FoV and for large base height ratio points. On the other hand, the "Fair" INS class exhibits larger WCD uncertainties (factor of 10) with higher variability in the obtained WCD uncertainty. The airborne observations generally result in larger confidence regions and higher WCD uncertainties compared to the drone scenario. An increased distance for the feature rays results in an increased dilution of the attitude uncertainty and therefore larger uncertainties are to be expected for higher flight altitudes.

Our results suggest that uncertainties in the airborne scenario are approximately six times greater than those in the drone scenario, despite a flight altitude ratio of 16:1 (2000 m for airborne versus 120 m for drone). Crucially, our study maintains consistent key parameters such as base-height ratio, FoV, and image overlap across both scenarios. In a purely geometric interpretation without probabilistic modeling, one might expect uncertainties to scale proportionally with flight altitude. In this regard, our findings indicate a statistical reduction in vertical uncertainty due to the intersection of two lines of sight, which is indicated in Figure 9 by the intersection likelihood compared to the single line-of-sight likelihood. We hypothesize that this effect increases with the altitude elucidating why increased distances in the airborne scenario do not linearly translate into increased uncertainty, highlighting the influence of probabilistic modeling in the variability of vertical uncertainties between different flight altitudes.

Regarding the viewing geometry influence, WCD uncertainties for both drone and plane scenarios follow a similar order (increasing from $P_r$ to $P_c$ to $P_l$) and show comparable patterns of variability.

Interestingly, the uncertainties for point $P_c$ (viewed by cross lines) exhibit high variability, especially in the airborne flight and the "Fair" INS case. In such a scenario, this increased variability of the uncertainty estimates for point $P_c$ can be attributed to the unique interplay between its viewing geometry and the anisotropic attitude noise with low concentration in the "Fair" INS case. The viewing geometry of $P_c$ inherently results in intersecting lines-of-sight but the sampled lines-of-sight, in the "Fair" INS class and high altitude airborne flight, results in greater variability in the positional likelihood near the MLE. This, coupled with the effect of refraction, can account for this pronounced variability in the WCD uncertainty estimates for $P_c$ viewing geometry.

4.2.2. Evaluation of Uncertainty Metrics

In this section, we present the results of our likelihood-based inference approach for estimating the WCD uncertainty using the $t_o$ and r test statistics. The rejection rates of the null hypothesis (True parameter value) were calculated through Monte Carlo simulations, employing a significance threshold of $\alpha =0.05$.

Figure 11 displays the rejection rates for various combinations of INS performances, depths, and viewing geometries in the drone scenarios. The statistical tests revealed a strong correlation among the samples, with rejection rates consistently aligning closely with the theoretical rejection rate of 5%, ranging from 3.5% to 7.0%. These findings indicate that both the $t_o$ and r test statistics are appropriate for estimating 95% CI for the WCD parameter. The observed correlation suggests that the WCD uncertainties are symmetric in the drone scenarios. Although the Expected Fisher information was not investigated in this study, the low variability of the uncertainties in most cases (orange error bars in Figure 10), indicates that it would have yielded similar results and performance to the Observed Fisher Information.

For the airborne scenarios, we observed a similar performance of the test statistics, except for point $P_c$, which showed distinct rejection rates, as illustrated in Figure 12.

Figure 11. Drone scenarios : rejection rates of the $t_o$ Wald statistic (Observed Fisher Information) and the signed likelihood ratio statistic r. The x-axis represents different simulated depths, while the y-axis represents the rejection rate percentage based on a CL of 95%. The green dashed horizontal line indicates the 5% rejection rate threshold.

Figure 11. Drone scenarios : rejection rates of the $t_o$ Wald statistic (Observed Fisher Information) and the signed likelihood ratio statistic r. The x-axis represents different simulated depths, while the y-axis represents the rejection rate percentage based on a CL of 95%. The green dashed horizontal line indicates the 5% rejection rate threshold.

Figure 12. Airborne scenarios : rejection rates of the $t_o$ Wald statistic (Observed Fisher Information) and the signed likelihood ratio statistic r. The x-axis represents different cases, while the y-axis represents the rejection rate percentage based on a confidence level of 95%. The green dashed horizontal line indicates the 5% rejection rate threshold.

Figure 12. Airborne scenarios : rejection rates of the $t_o$ Wald statistic (Observed Fisher Information) and the signed likelihood ratio statistic r. The x-axis represents different cases, while the y-axis represents the rejection rate percentage based on a confidence level of 95%. The green dashed horizontal line indicates the 5% rejection rate threshold.

Specifically, for the "Fair" INS scenario, the rejection rates for the $t_o$ test statistic exceeded 15% at point $P_c$, while the r statistic demonstrated relatively consistent rejection rates across different points (3.5%-7%). Given the WCD uncertainty results of $P_c$, highlighted in the previous section, we can infer that the higher variability in the WCD uncertainty estimates for $P_c$ in the previous subsection is consistent with the observed rejection rates. These findings suggest that the $t_o$-based CI may overstimate of WCD uncertainty in high altitude scenarios. This is particularly the case for viewing geometries similar to $P_c$ with a 19% rejection rate for the most noisy attitude scenario. On the other hand, the r statistic showed greater effectiveness in estimating the WCD uncertainties in this particular case ($P_c$ and airborne "Fair" INS).

In summary, the evaluation of uncertainty metrics through rejection rate analysis provides valuable insights into the reliability of uncertainty estimation methods for both the drone and airborne cases. Both scenarios exhibit reasonably consistent performance, except under extreme attitude noise and non-parallel lines viewing geometries.

4.3. Water Air Interface Height Inference

In our experimental setup, we sought to present the profile likelihood of h under multi-point and multi-view likelihood scenarios. To begin with, we need to clarify the conditions under which the WAI height is inferable from camera poses. Primarily, when the lines-of-sight are derived from parallel lines (i.e. $P_l$ or $P_r$ viewing geometry), or have equal incidence angles (i.e. $P_c$ viewing geometry), the lines-of-sight having the maximum likelihood will always intersect, irrespective of the WAI height. In these cases, the parameter h cannot be inferred because it has no influence on the likelihood. In contrast, when the lines-of-sight are non-coplanar with the baseline [28], then the lines-of-sight having the maximum likelihood will not intersect unless they are refracted under an optimal value of the WAI height parameter. This distinct viewing geometry can be readily obtained by introducing different incidence angles in the simulated lines-of-sight, allowing for the inference of the parameter h through the profile likelihood analysis. Either translating the same points ($P_c$, $P_l$ and $P_r$ ) or introducing different flight altitudes provide viewing geometries with different incidence angles. In our setup, points $P_l$, $P_r$ and $P_c$ are translated in both x and y dimensions, with different offsets ${\delta }_x$ and ${\delta }_y$ respectively. For drone simulations, ${\delta }_x=5m$ and ${\delta }_y=10m$, whereas for the airborne scenarios, these offsets are scaled by the altitude ratio. Applying these adjustments, we use the likelihood 20 for the WAI height parameter h inference with a WAI height true value of $\tilde{h}=0$, for varying sample sizes, namely $1,10,100,1000,10000$. Each sample refers to a couple of measured camera poses for a single feature point given the refraction at $\tilde{h}$. This means that for n samples, the parameter is of dimension $9*n+1$ (1 for h, $3n$ for ${\mathit{x}}_i$ and $6n$ for ${\mathit{x}}_{\mathit{s}}$).

Figure 13 shows on the left side WAI height normalized profile likelihood obtained for the translated point $P_c$, at a depth of 5 meters and for two sample sizes : 1 and 10. The associated r statistic is presented on the right side of the figure with the CI for two levels, namely 95% and 68%. According to the profile involving only one sample, (orange profile), two non-intersecting lines-of-sight can provide an optimal value for the WAI height parameter h. For both 1 sample and 10 samples profiles, we observe a typical asymmetry in the profile which is characterized by a flatness on the left side (below 5 meters depth). This pattern can be mainly attributed to the fact that when the WAI is relatively below the WCD for a given feature point, the positional likelihood should not decrease with the WAI height since there is no refraction above the WAI. The r statistic, affected by this flatness, successfully captures this asymmetry, which results in a non-identifiable 95% CI lower bound for the 1 sample inference although we could delineate 68% CI. On the other hand, 10 samples provides sufficient statistical information for delineating both 95% bounds for the WAI height h. However, based on the r 95% CI, the achieved precision of approximately 6.5 m for 10 samples is notably high. Therefore, it would be valuable to explore the possible improvement in WAI uncertainty when larger sample sizes are taken into account.

Figure 14 presents heatmaps of WAI height uncertainties, determined by r CI for point $P_c$ across different INS classes for the drone scenario. These heatmaps account for a variety of depths and sample sizes. As the sample sizes increased, so did the statistical information, leading to more precise profile likelihoods and subsequently, smaller uncertainties in the WAI height. Interestingly, a trend of uncertainty decreasing by a factor between 3 and $\sqrt{10}$ has been noted. This pattern, where uncertainty contracts in relation to the square root of the sample size ($\sqrt{n}$), echoes the square-root law common in Gaussian errors and is linked to the central limit theorem’s principle of sample averages approaching a normal distribution as the sample size expands, leading to decreased uncertainty. Analogously to the WCD inference, the "Good" and "Excellent" INS demonstrated similar performances. Intriguingly, the depth has no major impact on the WAI uncertainties except for the observed lag of the 1 meter depth scenario with sample size compared to other depths. This insinuates that the task of WAI inference tends to be more difficult when it relies on observations taken at depths near the WAI surface.

In Figure 15, we display the depth-averaged WAI uncertainties, together with their corresponding standard deviations, spanning viewing geometries, INS classes and sample sizes, for both drone and airborne settings. It’s often the case that the variability of WAI uncertainty with depth, as indicated by error bars, is more pronounced in difficult scenarios (10 samples). For the WAI height inference, the viewing geometry $P_c$ showed the best performance overall, suggesting that cross lines are strongly relevant in pushbroom acquisitions if one is interested in inferring the WAI height. It is particularly noteworthy that the achieved WAI height 95% CI of 17 cm and 19 cm for ’Excellent’ and ’Good’ quality INS, respectively, closely aligns with the Global Navigation Satellite System (GNSS) precision of 5 cm, when considering a substantial sample size of $10,000$. This suggests that with relatively small high resolution data sets, drone-based acquisitions can provide significantly accurate estimations of the WAI height in optically clear shallow waters with a textured seabed. Furthermore, considering the decline in uncertainties by a factor of $\sqrt{n}$ with increasing sample size, even with the "Fair" INS quality, we extrapolate that drone flights could achieve a WAI height precision of 16 cm given a sample size of about 1 million. Although 1 million matched feature points may seem substantial for a single survey of pushbroom imagery in coastal areas, given constraints like water clarity and a textured seabed, such a number of points is easily attainable with high-resolution frame video imagery under optimal conditions. This emphasizes the efficiency of drone-based acquisitions even with less-than-ideal INS quality, given an adequately large data set of paired feature points.

In contrast, airborne measurements yield a higher degree of uncertainty in WAI estimation with no successful inference for the "Fair" INS . Similarly to the WCD inference results, this is largely attributed to the increased operational altitude that intensifies optical distortion and viewing geometry variability due to a more distant and extensive footprint. Despite these challenges, airborne measurements with high quality INS achieve a metric precision for a 10,000 sample size, allowing to extrapolate a centimeter precision level for a 1 million sample size.

5. Discussion

The proposed study highlights the significant role of the camera pose quality on the WCD uncertainty. The "Excellent" INS class, typically deemed a high standard in hydrography, did not achieve satisfactory WCD precision (around 2 m) for airborne flights. This discrepancy underlines the advantages of drone-based surveys operating at lower altitudes with high quality camera pose measurements, which provided high WCD precisions, specifically around 50 cm for high base height ratios like $P_r$. According to International Hydrographic Organization (IHO) standards [29], this precision level falls within the ’Order 1a’ requirements without taking account of depth influence. However, the inherent limitations of higher base-height ratios and low flight altitudes, such as reduced and inefficient coverage, need to be considered. In contrast, for points $P_l$ and $P_c$, the uncertainty ranged between .7-1m, meeting the requirements of the ’Order 2’. Furthermore, we state that combining frame and pushbroom imagery could further improve camera pose quality for high altitude flights. Indeed, our geometric likelihood approach enables through-water photogrammetry bundle adjustment, presenting the opportunity to utilize frame imagery in order to improve camera pose accuracy. While our study provides key insights into camera pose uncertainties in through-water photogrammetry, it is important to note a limitation in our experimental design concerning multi-view geometry principles. Specifically, our experiments did not accommodate the principle of utilizing overlaps in both across and along flight directions, a common approach in multi-view geometry to augment camera pose samples and reduce uncertainty in triangulation. Our focus was on maintaining a uniform analysis across two camera poses in order to directly compare the effects of stereo geometry. This approach is also appropriate given that pushbroom studies often rely on stereo-pair configurations, where extensive cross-line datasets are very limited. As we focused on pushbroom geometric scenarios in our study, we only investigated two camera poses for triangulation, and we hypothesize that having an adequate number of feature points for each camera pose can refine the pose and eliminate the reliance on challenging-to-obtain GCP in coastal areas. Such advancements in camera pose quality have the potential to significantly enhance the precision levels of drone observations, potentially meeting stringent standards like the Special order in clear shallow waters with textured seabeds.

Furthermore, based on the WCD uncertainty results, the Observed Fisher Information performance and its low sensitivity to sampling suggest that classical variance-covariance propagation approaches, are reliable except for certain challenging scenarios. These include high flight altitudes, low quality camera pose measurements, and viewing geometries with non co-planar FoV ($P_c$ case). These viewing geometries are more prevalent in frame imagery since the spatial footprint in frame imagery is rectangular, thus establishing the likelihood ratio statistic as a potential tool for robust uncertainty estimation in such scenarios.

Beyond the evaluated uncertainties, it’s worth emphasizing that our results contribute to a deeper understanding of photogrammetric techniques for coastal photogrammetric surveying. The ability to accurately infer the water surface position from imagery alone reduces the reliance on additional tide data or separate WAI height observations. Although our initial assumptions considered a horizontal water surface, not accounting for surface fluctuations, the results showed satisfactory precision of around 17 cm when using large samples and non-intersecting viewing geometries. This precision level is particularly impressive considering calm sea conditions and light weather situations, which can be encountered in sheltered area with low significant wave heights. Nonetheless, it’s important to consider that in a real-world scenario, coastal waters fluctuations could occur due to a range of factors such as changes in tide and weather, which may influence the temporal coherence of the imagery dataset. This consideration indicates that our approach can be further improved by taking into account these potential fluctuations. For instance, simulations of wave spectra and sea state conditions can be used to create more realistic scenarios, offering a comprehensive perspective on the influence of surface fluctuations on WAI height uncertainties and their subsequent impact on the derived bathymetric estimations. The noted viewing geometry condition for WAI height inference might favor in-track stereo pairs more then across-track satellite stereo pairs which are often acquired with the same incidence angles. Therefore, a continuation of this study would be extending the investigation of the effect of viewing geometry and proposing appropriate flight acquisition modes for WAI height inference. Also, the performance of r confidence intervals that we utilized to estimate uncertainties in WAI height were not assessed through hypothesis testing, in contrast to what we performed for the WCD. This examination was out of the scope of this study and therefore, could serve as a potential area for further investigation and research. Moreover, our approach for WAI inference can be effectively adapted with stationary high-quality stereo camera observations as a technique for tide monitoring.

The results also invite a deeper examination of innovative uncertainty modeling methods. Such methods can be beneficial in harmonizing the geometric inference of WCD with other forms of inference like radiometric estimation with analytical radiative transfer models, within a unified theoretical likelihood framework as explored by Jay et al. [30], Sicot et al. [31]. Interestingly, the relative insensitivity of the geometric inference to depth variations, which contrasts with its radiometric counterpart as suggested by these studies, indicates potential complementarity between the geometric and radiometric inference. However, it’s important to recognize that both approaches, geometric and radiometric, require a degree of optical water clarity. Therefore, the impact of depth can introduce inherent constraints on both types of inferences, especially in environments with turbid waters.

Our study is a robust step forward in uncertainty modeling for bathymetric photogrammetry in coastal shallow waters and has unveiled several insights. The main focus of our study is the investigation of camera pose uncertainties in through-water photogrammetry. However, it is important to acknowledge that uncertainties may arise from additional sources, such as matching errors, camera imperfections and wave dynamics, which are significant factors in coastal waters imagery.

6. Conclusions

This study represents a detailed exploration of the potential and accuracy of through-water photogrammetry, with a specific emphasis on the WCD and the WAI height inferences. Our observations reveal the profound influence of viewing geometry and camera pose quality on the resulting uncertainties, overshadowing the impact of depth.

Our approach bridges the gap between advanced probabilistic modeling and stereo-photogrammetric triangulation. This innovative integration provides a comprehensive understanding of the complexities associated with through-water 3D reconstruction methods. The utility of drone technology, complemented by high quality camera pose measurement and spectral imagery, proves to be a compelling tool for high-precision through-water photogrammetry. Importantly, these advancements contribute to the provision of a framework for stereo-photogrammetry predictive uncertainties in bathymetric charting. Additionally, we demonstrate that inferring the water surface elevation is achievable from camera pose measurements reducing the dependence on auxiliary data sources such as tide level measurements or WAI height observations.

Notwithstanding the success of our investigations, we acknowledge the limitations of our current study. Specifically, the questions related to providing compelling recommendations for viewing geometries and hypothesis testing for WAI height, requiring further exploration and elucidation. To enhance the robustness and reliability of our methodology under a variety of conditions, we encourage future research to incorporate surface fluctuations and camera optical imperfections.

In conclusion, our study sets a robust foundation for future investigations in the field of through-water photogrammetry. It contributes significantly towards refining photogrammetric techniques, and is shaping the future of through-water surveying and data collection methodologies in coastal areas. Although this study marks only the initial steps in this evolving field, we are confident that our efforts have established an essential milestone, paving the way for further advancements in through-water photogrammetry.