Revolutionizing Submersible Safety through Advanced Mathematical Modeling: From Bayesian Searching Model to Dynamic Oceanic Simulations

1. Introduction

Submersibles offer unique opportunities for underwater exploration, but their safe operation faces challenges such as loss of communication, mechanical failures, and dynamic marine environments. Accurate prediction of submersible locations and effective emergency procedures are crucial for ensuring the safety of both passengers and crew.

This paper aims to:

• Develop a Seawater Density Model to represent depth-dependent density variations.
• Construct a Submersible Mechanical Model for trajectory prediction.
• Propose a Bayesian Searching Model and Extended Kalman Filter for locating missing submersibles.
• Perform sensitivity analysis and extrapolate models to other environments.

2. Locate

To precisely locate the submersible, we developed two key models: the Seawater Density Model and the Submersible Mechanical Model. These models allow us to accurately determine the submersible’s position after losing communication by analyzing both the surrounding environment and the dynamics of the submersible itself.

2.1. Seawater Density Model

The first step in predicting the movement of the submersible is to analyze the environment, specifically the seawater density, which determines both the buoyancy and neutral buoyancy depth of the submersible. This is essential for understanding the submersible’s equilibrium behavior.

2.1.1. Seawater Density Stratification

The density of seawater varies due to thermohaline stratification, which refers to the changes in temperature and salinity at different depths [3]. Surface waters are generally warmer and less saline, leading to lower density, while deep waters are colder, saltier, and therefore denser.

Based on literature [2] and our assumptions, we model seawater density as a hyperbolic tangent function of depth, $\rho \left(z\right)$. The density generally increases with depth due to higher pressure and changes in temperature and salinity, resulting in a nonlinear relationship.

$$\rho \left(z\right)={\rho }_0+\sum _{i=1}^n{A}_i·\left(tanh\left({\displaystyle \frac{z-z_{0i}}{B_i}}\right)-1\right)$$

(1)

Here, ${\rho }_0$ represents the sea-level water density, $2A_i$ represents the density increment for each layer, and $B_i$ controls the steepness of the density gradient.

Figure 1. Fitted density stratification of seawater.

Figure 1. Fitted density stratification of seawater.

Using this model, we stratify seawater into three layers with distinct pycnoclines. The mean density of seawater in our study area is found to be 1039.67 kg/m³.

2.2. Submersible Mechanical Model

After analyzing the seawater density, we then model the submersible itself. This involves both static and dynamic analysis to understand how the submersible moves in response to various forces.

2.2.1. Static Buoyancy Analysis

We first analyze the equilibrium state of the submersible. Let the submersible’s initial position and velocity at the time of losing communication be $(x_0,y_0,z_0)$ and $(v_r,v_z)$, respectively. We assume the mass of the submersible is constant (m) and that the propulsion force (F) has an angle $\gamma$ with respect to the x-y plane. The static buoyancy equation is given by:

$$Fsin\gamma +mg=mg·{\displaystyle \frac{\rho \left(z\right)}{{\rho }_s}}$$

(2)

Here, ${\rho }_s$ is the density of the submersible, and $\rho \left(z\right)$ is the density of the seawater. Solving this equation provides the neutral buoyancy depth $Z=Z_0$.

2.2.2. Dynamic Analysis

We then perform a dynamic analysis of the submersible, considering its movement and environmental factors. Figure 2 illustrates the free body diagram and velocity decomposition.

The submersible is subject to several forces, including propulsion ($\overrightarrow{F}$), gravity ($\overrightarrow{G}$), buoyancy (${\overrightarrow{F}}_{buo}$), and fluid resistance ($\overrightarrow{F_D}$):

$$m{\displaystyle \frac{d\overrightarrow{v}}{dt}}=\overrightarrow{F}+\overrightarrow{G}+{\overrightarrow{F}}_{buo}+{\overrightarrow{F}}_D$$

(3)

We decompose the velocity into components, where $\alpha$ is the angle between the velocity projection on the x-y plane and the x-axis, and $\beta$ is the angle between the velocity and the x-y plane. $\gamma$ and $\delta$ represent the angles of the propulsion force with respect to the x-y plane and the x-axis.

The system of differential equations representing the submersible’s dynamics is given by:

$$\left\{\begin{array}{c}Fsin\gamma -{\displaystyle \frac{1}{2}}\rho \left(z\right)v^2{C}_{D}Asin\beta -{\displaystyle \frac{\left(\rho \left(z\right)-{\rho }_0\right)}{{\rho }_0}}mg=m{\displaystyle \frac{d{v}_z}{dt}}\hfill \\ Fcos\gamma cos\delta -{\displaystyle \frac{1}{2}}\rho \left(z\right)v^2{C}_{D}Acos\beta cos\alpha =m{\displaystyle \frac{d{v}_x}{dt}}\hfill \\ Fcos\gamma sin\delta -{\displaystyle \frac{1}{2}}\rho \left(z\right)v^2{C}_{D}Acos\beta sin\alpha =m{\displaystyle \frac{d{v}_y}{dt}}\hfill \end{array}\right.$$

(4)

Table 1. Dynamic parameters for different submersibles.

Parameter	Submersible 1	Submersible 2
Velocity $v_0$ (m/s)	15	30
Propulsion force $F_0$ (N)	28	34
Mass m (kg)	1000	1000
Drag coefficient $C_D$	0.5	0.5
Reference area A (${\mathrm{m}}^2$)	0.6	0.6
${\rho }_s$ (${\mathrm{kg}/\mathrm{m}}^3$)	1022.75	1030.75

Table 1. Dynamic parameters for different submersibles.

Parameter	Submersible 1	Submersible 2
Velocity $v_0$ (m/s)	15	30
Propulsion force $F_0$ (N)	28	34
Mass m (kg)	1000	1000
Drag coefficient $C_D$	0.5	0.5
Reference area A (${\mathrm{m}}^2$)	0.6	0.6
${\rho }_s$ (${\mathrm{kg}/\mathrm{m}}^3$)	1022.75	1030.75

2.2.3. Boundary Conditions: Seabed Depth

The submersible’s motion is constrained by the seabed, which limits the depth the submersible can reach. Figure 3 shows the seabed depth data for the Ionian Sea [1,5].

We solved the system of differential equations using the 4th order Runge-Kutta method to accurately predict the trajectory of the submersible after losing communication. Figure 4 and Figure 5 depict two scenarios: one where the submersible hits the seabed and one where it maintains neutral buoyancy.

2.3. Uncertainties in the Model

2.3.1. Uncertainty in Mechanical Default

It is challenging to determine whether mechanical failure occurs before or after a communication loss. A monitoring system for key components can help detect failures early and trigger alarms or automated measures.

2.3.2. Noise Uncertainties and Solutions

Noise from fluid dynamics, such as ocean currents, affects the submersible’s trajectory. To improve accuracy, GPS and IMU sensors are used to provide position $(x,y,z)$ and velocity $(v_x,v_y,v_z)$ data, respectively. These sensors, however, introduce noise. To reduce uncertainty, we apply the Extended Kalman Filter (EKF), which is suitable for the non-linear nature of submersible motion.

The EKF prediction and update steps are [4]:

$${\widehat{\mathit{x}}}_{k|k-1}=f({\widehat{\mathit{x}}}_{k-1|k-1},{\mathit{u}}_{k-1}),{\widehat{\mathit{x}}}_{k|k}={\widehat{\mathit{x}}}_{k|k-1}+{\mathit{K}}_k{\tilde{\mathit{y}}}_k$$

(5)

where ${\mathit{K}}_k$ is the Kalman gain and ${\tilde{\mathit{y}}}_k$ is the measurement innovation.

The state transition equations for submersible motion are:

$$\begin{array}{cc}\hfill v_{x,t|t}& =v_{x,t-1}+{\displaystyle \frac{\Delta t}{m}}\left(Fcos\gamma cos\delta -{\displaystyle \frac{1}{2}}\rho \left(z\right)C_{D}A{v}_x\right)\hfill \\ \hfill v_{y,t|t}& =v_{y,t-1}+{\displaystyle \frac{\Delta t}{m}}\left(Fcos\gamma sin\delta -{\displaystyle \frac{1}{2}}\rho \left(z\right)C_{D}A{v}_y\right)\hfill \\ \hfill v_{z,t|t}& =v_{z,t-1}+{\displaystyle \frac{\Delta t}{m}}\left(Fsin\gamma -{\displaystyle \frac{1}{2}}\rho \left(z\right)C_{D}A{v}_z-{\displaystyle \frac{\rho \left(z\right)-{\rho }_0}{{\rho }_0}}mg\right)\hfill \end{array}$$

(6)

2.3.3. Enhancement Model: BIKF and DKF Gaussian Noise Mixture Model

To simulate unpredictable ocean currents, we use a Gaussian Mixture Model (GMM) for non-Gaussian noise, combining small- and large-variance Gaussian noise. The Background Impulse Kalman Filter (BIKF) handles non-Gaussian noise more effectively [4]. We also apply a Distributed Kalman Filter (DKF) for data fusion [6].

3. Searching

After locating the location of the submersible, we need create a model utilizing data from previous models used in location task to suggest starting points and searching trajectories for the equipment, aiming to reduce the time needed to locate a missing submersible. We need to determine the likelihood of locating the submersible based on time and the accumulated search results.

Here, we propose a Bayesian Searching Model to solve the problem.

3.1. Bayesian Searching Model

We can have a prior probability about the submersible’s location, which can be obtained from the positioning in the first section. However, each time the rescue boat conducts a search, it brings new information about the location of the submersible. To fully utilize this information, we employ a Bayesian search model. We use the results from each search as the prior probability distribution for our next search[?]. This iterative process is aimed at minimizing the search time as much as possible.

3.1.1. Bimodal Gaussian Distribution Function

Due to the model uncertainties mentioned in the first section, we cannot determine whether the mechanical defects will happen after losing communication. So we need to consider possibilities in both situations into consideration, mechanical defects happens and no mechanical defects after signal lost. We use bimodal Gaussian Distribution to incorporate the these two situations, resembling the prior probability density of searching.

Suppose the last signal before losing communication transmit the state of submersible is $(x,y,z,v_x,v_y,v_z)$. The state include the position and velocity in different directions.

From the previous analysis, the submarine will quickly enter a static state in the z-direction. Consider the influence of ocean currents and noise. we consider two situations where the submersible continues to sail and mechanical failure, representing the prior probability map at the initial moment as a bimodal Gaussian distribution [8], that is

$$f(x,y)=w_1·{\displaystyle \frac{1}{2\pi {\sigma }_1^2}}exp\left(-{\displaystyle \frac{{(x-x_1)}^2+{(y-y_1)}^2}{2{\sigma }_1^2}}\right)+w_2·{\displaystyle \frac{1}{2\pi {\sigma }_2^2}}exp\left(-{\displaystyle \frac{{(x-x_2)}^2+{(y-y_2)}^2}{2{\sigma }_2^2}}\right),$$

(7)

${\omega }_1$ and ${\omega }_2$ are respectilythe possibilities of mechanical defects and the no mechanical defects. ${\sigma }_1$ and ${\sigma }_2$ are the influence of noise and current in two situations respectively.

Table 2. Parameters in bimodal Gaussian distribution function.

parameters	value	parameters	value
${\omega }_1$	0.7	${\omega }_2$	0.3
${\sigma }_1$	10	${\sigma }_2$	10
$x_1$	3800	$y_1$	3800
$x_2$	2500	$y_2$	2500

Table 2. Parameters in bimodal Gaussian distribution function.

parameters	value	parameters	value
${\omega }_1$	0.7	${\omega }_2$	0.3
${\sigma }_1$	10	${\sigma }_2$	10
$x_1$	3800	$y_1$	3800
$x_2$	2500	$y_2$	2500

Using the bimodal Gaussian distribution function in eq (7), we select sea area of $6\mathrm{km}\times 6\mathrm{km}$, the probability distribution graph of submersible’s location (2-dimensional) is shown below.

Figure 6. The probability distribution graph of submersible location.

Figure 6. The probability distribution graph of submersible location.

3.1.2. Discretization of Searching Possibilities (Avocado Grid Searching)

We assume the effective detection range of the rescue boat conducted is within time interval T. Assume in time interview T, the unit searching area of rescue ship is D. So we can divide the sea into grids with an area of D in our model. Which looks just like an avocado.

These grids are regarded as n non-interfering areas $R_i$($i=1,2,...,n$). Use $x_{i,k}$ to represent the status in the ith area at the kth moment. The probability of each grid is the integral of the probability density in the area of the grid. $y_{i-}$, $y_{i+}$, $x_{i-}$, $x_{i+}$ are the boundaries of the grid.

$$p(i,t)={\int }_{y_{i-}}^{y_{i+}}{\int }_{x_{i-}}^{x_{i+}}f(x,y)dxdy$$

(8)

After calculating the probability for each grid through integration, we normalize these values for each grid point.

The probability of the sensor successfully observing the submersible at its location is denoted as q (conditional probability). p is the probability of submersible located in current observed grid. Based on Bayes’ theorem, we derived the following formula for updating the posterior probability with each combined observation: After this detection, if the submersible is not detected in this grid, the possibility of this grid under detected is renewed.

Figure 7. Possibilities renew after current detection.

Figure 7. Possibilities renew after current detection.

After detection, the possibility of the grid being detected is renew to $p^{\prime }$

$$p^{\prime }={\displaystyle \frac{p(1-q)}{(1-p)+p(1-q)}}={\displaystyle \frac{1-q}{1-pq}}p{r}^{\prime }={\displaystyle \frac{1}{1-pq}}r$$

(9)

After detection, the possibility of the rest of the grids is $r^{\prime }$ :

$$r^{\prime }={\displaystyle \frac{1}{1-pq}}r$$

(10)

Using eq (9,10), we continue to iterate the searching process, combining with the submersible mechanical model in the first section, we have the simulation of searching as Figure 8. In our case, the bimodal Gaussian distribution has a large variance. When the variance is small, a bimodal image will appear. When the variance is large, it will produce the image of an avocado that we speak of.

The previous graph shows the iteration of posterior possibilities of submersible’s location. Combining the model in locating and the searching we simulate the real search route based on the dynamic equations and Bayesian searching model of the submersible. These graphs show the MSE between different predicted trajectories and real trajectory.

Current we are simulating the overlapping process of the movement of submersible and the searching route. The red line in Figure 9(a) simulate the real movement of submersible on x-y plane. In the process of Bayesian searching iteration and the rescue boat searching, when the route of submersible is overlapped with currently searching grid, the successful searching possibilities is q, which is the largest.

3.2. The Probability of Finding the Submersible as a Function of Time and Accumulated Search Results

The probability of locating the submarine is expressed as a function of time and accumulated search results. Let $H_i$ be a binary variable used to determine whether the submarine is found in the ith search.

$$H_i=\left\{\begin{array}{cc}1\hfill & \mathrm{if}\mathrm{found},\hfill \\ 0\hfill & \mathrm{if}\mathrm{not}\mathrm{found}\hfill \end{array}\right.$$

(11)

Every searching trial has two results, the submersible either being searched or not being searched. Suppose at the nth detection, the probability of the submersibles located in this searching grid is $p_n$. where $p_n$ is derived from the prior probability before each detection.

In time $i\ast t$, the probability of submersible being successfully searched is

$$P\left(\sum _{n=1}^i{H}_i\right)=1$$

(12)

Notice the ∑ summation here means or. Moreover, because of

$$P(\sum _{n=1}^i{H}_i=1)=1-P(\sum _{n=1}^i{H}_i=0)=1-\prod _{n=1}^{i}P(H_n=0)=1-\prod _{n=1}^i(1-p_n·q)$$

(13)

In summary, the probability of finding the submersible as a function of time and accumulated search results can be expressed as:

$$P(\sum _{n=1}^i{H}_i=1)=1-\prod _{n=1}^i(1-p_n·q)$$

(14)

Here the q is set to be 0.9. The probability of finding the submersible over searching trials will finally approach q.

Table 3. Probability of finding the submersible.

Searching trial	20	40	60	80	100
P(${\sum }_{n=1}^i{H}_i=1$)	0.4065	0.6352	0.7611	0.8309	0.8763

Table 3. Probability of finding the submersible.

Searching trial	20	40	60	80	100
P(${\sum }_{n=1}^i{H}_i=1$)	0.4065	0.6352	0.7611	0.8309	0.8763

4. Extrapolation

The models presented earlier need to be adapted to account for varying environmental conditions and extended for the detection of multiple submersibles in close proximity. As an example, we consider the Caribbean Sea.

4.1. Seawater Density in the Caribbean Sea

According to oceanographic data, the seawater density in the Caribbean Sea remains relatively consistent, ranging from 1023.5 to 1024 $\mathrm{kg}/{\mathrm{m}}^3$ [?]. In contrast to previous models where the $tanh$ function was employed to represent seawater density, the $tanh$ function is not suitable for the Caribbean Sea. Therefore, we approximate seawater density using a linear function based on the minimum riverbed elevation:

$$\rho \left(z\right)=\left(\rho \left(z_{\max }\right)-\rho \left(0\right)\right)·{\displaystyle \frac{z}{z_{\max }}}+\rho \left(0\right)$$

(15)

where $\rho \left(z_{\max }\right)$ is the maximum seawater density, and $\rho \left(0\right)$ is the seawater density at the surface.

Figure 10. Probability of finding the submersible over searching trials.

Figure 10. Probability of finding the submersible over searching trials.

4.2. Warning and Obstacle Avoidance System for Multiple Submersibles

To enable multiple submersibles to operate safely within the same maritime area, we have implemented a warning and obstacle avoidance system. The warning distance between two submersibles is denoted as $d_w$.

Figure 11. Schematic diagram of obstacle avoidance principle.

Figure 11. Schematic diagram of obstacle avoidance principle.

Considering two submersibles with trajectories ${\overrightarrow{r}}_1\left(t\right)$ and ${\overrightarrow{r}}_2\left(t\right)$, respectively, the submersibles will adjust their propulsion forces whenever their predicted trajectories approach a dangerous proximity, i.e., $\left|\right|{\overrightarrow{r}}_1\left(t\right)-{\overrightarrow{r}}_2\left(t\right)\left|\right|<d_w$. The submersibles adjust their paths while returning to their neutral buoyancy levels to avoid collision:

$$\left\{\begin{array}{c}{\overrightarrow{F}}_1(t+1)={\displaystyle \frac{C_1·F_w}{\left|\right|{\overrightarrow{r}}_1\left(t\right)-{\overrightarrow{r}}_2\left(t\right)\left|\right|}}·{\displaystyle \frac{({\overrightarrow{r}}_1\left(t\right)-{\overrightarrow{r}}_2\left(t\right))}{\left|\right|{\overrightarrow{r}}_1\left(t\right)-{\overrightarrow{r}}_2\left(t\right)\left|\right|}}\hfill \\ {\overrightarrow{F}}_2(t+1)={\displaystyle \frac{C_2·F_w}{\left|\right|{\overrightarrow{r}}_2\left(t\right)-{\overrightarrow{r}}_1\left(t\right)\left|\right|}}·{\displaystyle \frac{({\overrightarrow{r}}_2\left(t\right)-{\overrightarrow{r}}_1\left(t\right))}{\left|\right|{\overrightarrow{r}}_2\left(t\right)-{\overrightarrow{r}}_1\left(t\right)\left|\right|}}\hfill \end{array}\right.$$

(16)

Figure 12. Comparisons of submersible trajectories in the Caribbean Sea.

Figure 12. Comparisons of submersible trajectories in the Caribbean Sea.

Figure 13(a) illustrates the trajectories of submersibles with different densities under external forces, while Figure 13(b) depicts the depth variation over time. The densities of the three submersibles are 1.38, 1.40, and 1.42, respectively.

The highlighted regions indicate moments when the submersibles were influenced by sudden external forces, such as rapid ocean currents, which caused a temporary deviation from their neutral buoyancy levels. Despite these disturbances, the submersibles quickly returned to their original buoyancy level, demonstrating the robustness of the model. Thus, the model exhibits stability in response to changes in density and short-term external forces, both over short and long time periods.

5. Reliability Analysis

In evaluating the reliability of the model, it is important to consider the assumptions made to simplify the model while ensuring it remains robust enough to provide useful insights in realistic scenarios. The assumptions provide a foundation for the model’s structure and allow for effective simulation under various conditions. Below is an analysis of the reliability of the model based on the key assumptions:

5.1. Assumptions and Their Impact on Reliability

• The shape and density of the submersible remain unchanged underwater.
  Impact on Reliability: This assumption simplifies the model by treating the submersible as a rigid object, unaffected by external pressures. In reality, while extreme water pressure could potentially cause minor deformations, the materials used in submersibles are generally designed to withstand these conditions. By ignoring minor shape deformations, we prioritize model simplicity without sacrificing significant accuracy, making this assumption both reasonable and supportive of model reliability.
• The density of seawater is influenced only by underwater depth.
  Impact on Reliability: By considering seawater density as dependent solely on depth, we eliminate the need to account for small changes in density caused by variations in longitude and latitude. This assumption holds because vertical movement significantly impacts the submersible’s buoyancy, while lateral movement over short distances has a negligible effect. This allows the model to focus on the primary variable—depth—thus increasing computational efficiency without compromising reliability for short timeframes and localized areas.
• Two distinct scenarios are assumed after losing communication: mechanical defect or no defect.
  Impact on Reliability: Simplifying the analysis to two scenarios—either the submersible continues with unchanged propulsion or propulsion ceases due to a mechanical defect—provides a straightforward framework for determining the submersible’s trajectory after losing communication. While other failure modes could exist in reality, this binary distinction allows the model to focus on the most critical outcomes, enhancing reliability in failure prediction without introducing unnecessary complexity.

5.2. Analysis of Model Reliability Based on Assumptions

• Structural Stability of the Submersible:
  The assumption that the submersible’s shape and density remain constant is integral to predicting reliable movement patterns underwater. Ignoring minor deformations caused by water pressure ensures that the model does not get bogged down by unnecessary complexities while maintaining realistic trajectory predictions. Given the robust materials used in submersibles, this assumption does not significantly reduce the model’s reliability.
• Seawater Density and Depth Considerations:
  The assumption that seawater density is only affected by depth simplifies the buoyancy calculation, which is crucial for predicting vertical movement. Because vertical movement has a more pronounced effect on submersible buoyancy compared to lateral movement, the model remains highly reliable in predicting submersible depth changes. While minor fluctuations in density based on geographic location are disregarded, this does not diminish the model’s ability to operate in small to medium-scale geographic areas over short periods.
• Failure Scenario Assumptions:
  The assumption of only two post-communication loss scenarios—continued propulsion or mechanical failure—enables a clear focus on the most likely outcomes. This binary framework is well-suited for simulating search and rescue operations. However, the model could be enhanced by incorporating more detailed failure modes, such as partial propulsion or sensor malfunctions. Nonetheless, the assumption allows for an efficient initial response framework and enhances the model’s reliability in the early stages of submersible tracking.

5.3. Overall Model Reliability

Based on these assumptions, the model remains reliable under typical operational conditions, particularly when used for short-duration missions or localized tracking of submersibles. By focusing on primary variables such as depth, buoyancy, and mechanical failures, the model provides robust predictions while maintaining computational simplicity.

Additional reliability could be gained by refining the failure scenarios or incorporating more real-world

6. Conclusions

In this paper, we presented a comprehensive mathematical model for enhancing submersible safety, focusing on location prediction, emergency procedures, and operational robustness. We developed and integrated models, including the Seawater Density Model, Submersible Mechanical Model, Bayesian Searching Model, and Extended Kalman Filter, to simulate and analyze the movement of submersibles in marine environments. These models were extended to different scenarios, such as the Caribbean Sea, and were equipped with a warning and obstacle avoidance system to manage multiple submersibles.

The reliability of the model was examined through sensitivity analysis, demonstrating its adaptability to varying environmental conditions and sudden external forces. The model’s accuracy was further improved through the incorporation of uncertainties using an Extended Kalman Filter, Gaussian Mixture Model, and Background Impulse Kalman Filter, ensuring precise trajectory prediction.

Incorporating both theoretical and practical elements, the model provided efficient search procedures through Bayesian search techniques, minimizing the time needed to locate a missing submersible. The study highlights the significance of simplifying assumptions while maintaining reliability, and the overall robustness of the models ensures their applicability across a range of marine environments and operational challenges.

Future improvements can focus on refining the failure scenarios, increasing the precision of real-time data inputs, and expanding the model to accommodate more complex oceanographic and mechanical dynamics.