1. Introduction

Table 1. Parameter Symbols and Descriptions.

2. Related Works

3. Methodology

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.

4. Results

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.

5. Discussion

6. Conclusions

Abbreviations

References

1. Singh, G.; Cottrell, R.; Eddy, T.; Cisneros-Montemayor, A. Governing the Land-Sea Interface to Achieve Sustainable Coastal Development. Frontiers in Marine Science 2021. [CrossRef]
2. Martínez, M.; Intralawan, A.; Vázquez, G.; Pérez-Maqueo, O.; Sutton, P.; Landgrave, R. The coasts of our world: Ecological, economic and social importance. Ecological Economics 2007, 63, 254–272. Ecological Economics of Coastal Disasters, . [CrossRef]
3. Pasquale, D.; Giovanni, B.R.; Vincenzo, P. Measurement for the sea; Springer Series in Measurement Science and Technology Series Editors, 2022.
4. Lyzenga, D. Passive remote sensing techniques for mapping water depth and bottom features. Applied Optics 1978, 379. [CrossRef]
5. Stumpf, R.; Holderied, K.; Sinclair, M. Determination of water depth with high-resolution satellite imagery over variable bottom types. Limnolology and Oceanography 2003, 48, 547,557. [CrossRef]
6. Lee, Z.; Carder., K.L.; Mobley, C.D.; Steward, R.G.; Patch, J.S. Hyperspectral remote sensing for shallow waters: 2. Deriving bottom depths and water properties by optimization. Applied Optics 1999, 38, 3831–3843. [CrossRef]
7. Hedley, J.; Roelfsemab, C.; Koetzc, B.; Phinn, S. Capability of the Sentinel 2 mission for tropical coral reef mapping and coral bleaching detection. Remote Sensing of Environment 2012. [CrossRef]
8. Cao, B.; Fang, Y.; Jiang, Z.; Gao, L.; Hu, H. Shallow water bathymetry from WorldView-2 stereo imagery using two-media photogrammetry. European Journal of Remote Sensing 2019, 52, 506–521. [CrossRef]
9. Mandlburger, G. Through-water dense image matching for shallow water bathymetry. Photogrammetric Engineering and Remote Sensing 2019, 85, 445–454. [CrossRef]
10. Hodúl, M.; Bird, S.; Knudby, A.; Chénier, R. Satellite derived photogrammetric bathymetry. ISPRS Journal of Photogrammetry ,Remote Sensing 2018, 142, 268–277. [CrossRef]
11. Murase, T.; Tanaka, M.; Tani, T.; Miyashita, Y.; Ohkawa, N.; Ishiguro, S.; Suzuki, Y.; Kayanne, H.; Yamano, H. A Photogrammetric Correction Procedure for Light Refraction Effects at a Two-Medium Boundary. Photogrammetric Engineering and Remote Sensing 2008, 74, 1129–1136. [CrossRef]
12. Agrafiotis, P.; Georgopoulos, A. Camera Constant in the Case of Two Media Photogrammetry. The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences 2015, XL-5/W5, 1–6. [CrossRef]
13. JCGM. Evaluation of measurement data - Guide to the expression of uncertainty in measurement, 2008.
14. Dolloff, J.T. RPC uncertainty parameters: Generation, application, and effects. Proceedings of the ASPRS Annual Convention, 2012, pp. 19–23.
15. Dolloff, J.; Theiss, H. The Specification and Validation of Predicted Accuracy Capabilities for Commercial Satellite Imagery. Proceedings of the ASPRS Annual Convention, 2014.
16. Slocum, R.K.; Wright, W.; Parrish, C.; Costa, B.; Sharr, M.; Battista, T.A. Guidelines for Bathymetric Mapping and Orthoisubsectionmage Generation using sUAS and SfM, An Approach for Conducting Nearshore Coastal Mapping. NOAA technical memorandum NOS NCCOS 265, United States, National Ocean Service; National Centers for Coastal Ocean Science (U.S.); United States, National Oceanic and Atmospheric Administration; Coral Reef Conservation Program (U.S.), 2019. [CrossRef]
17. Zhang, C.; Sun, A.; Hassan, M.A.; Qin, C. Assessing Through-Water Structure-from-Motion Photogrammetry in Gravel-Bed Rivers under Controlled Conditions. Remote Sensing 2022, 14. [CrossRef]
18. Sims-Waterhouse, D.; Isa, M.; Piano, S.; Leach, R. Uncertainty model for a traceable stereo-photogrammetry system. Precision Engineering 2020, 63, 1–9. [CrossRef]
19. Kalacska, M.; Lucanus, O.; Sousa, L.; Vieira, T.; Arroyo-Mora, J.P. Freshwater Fish Habitat Complexity Mapping Using Above and Underwater Structure-From-Motion Photogrammetry. Remote Sensing 2018, 10. [CrossRef]
20. Kiraci, A.C.; Toz, G. Theoretical analysis of positional uncertainty in direct georeferencing. The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences 2016, XLI-B1, 1221–1226. [CrossRef]
21. Skarlatos, D.; Agrafiotis, P. A novel iterative water refraction correction algorithm for use in structure from motion photogrammetric pipeline. Journal of Marine Science and Engineering 2018, 6. [CrossRef]
22. Tewinkel, G.C. Water Depths from Aerial Photographs. Photogrammetric Engineering 1963, pp. 1038–1042. Paper presented at the 29th Annual Meeting of the American Society of Photogrammetry, Washington, D. C.
23. Rofallski, R.; Luhmann, T. An Efficient Solution to Ray Tracing Problems in Multimedia Photogrammetry for Flat Refractive Interfaces. Journal of Photogrammetry Remote Sensing and Geoinformation Science 2022, 90. [CrossRef]
24. Mulsow, C. A Flexible Multi-media Bundle Approach. International Archives of Photogrammetry, Remote Sensing and Spatial Information Sciences, 2010, Vol. XXXVIII, pp. 472–477.
25. Nasiri, S.M.; Hosseini, R.; Moradi, H. The optimal triangulation method is not really optimal. IET Image Processing 2023, 17, 2855–2865. [CrossRef]
26. Bingham, C. An Antipodally Symmetric Distribution on the Sphere. The Annals of Statistics 1974, 2, 1201–1225.
27. Brazzale, A.; Davison, A.; Reid, N. Applied Asymptotics: Case Studies in Small-Sample Statistics; Cambridge Series in Statistical and Probabilistic Mathematics, Cambridge University Press, 2007.
28. Rinner, K. Problems of Two-Medium Photogrammetry. Photogrammetric Engineering 1948, pp. 275–282. Presented at the Annual Convention of the American Society of Photogrammetry, Washington, D. C.
29. International Hydrographic Organization. International Hydrographic Organization Standards for Hydrographic Surveys S-44 Edition 6.0.0, 2020.
30. Jay, S.; Guillaume, M.; Chami, M.; Audrey Minghelli, Y.D.; Lafrance, B.; Serfaty, V. Predicting minimum uncertainties in the inversion of ocean color geophysical parameters based on Cramer-Rao bounds. OPTICS EXPRESS 2018, 26. [CrossRef]
31. Sicot, G.; Ghannami, M.A.; Lennon, M.; Loyer, S.; Thomas, N. Likelihood Ratio statistic for inferring the uncertainty of satellite derived bathymetry. IEEE WHISPER 2021. [CrossRef]
32. Kay, L.S.M. Fundamentals of Statistical Signal Processing: Estimation Theory; Prentice Hall Signal Processing Series, 1993.
33. Severini, T.A. Likelihood Methods in Statistics; Oxford Statistical Series, 2000.