Method for predicting depth-averaged current velocities of underwater gliders based on data feature analysis

Underwater gliders are a new type of autonomous underwater vehicle (AUV). Compared with the traditional AUVs, underwater gliders are of low cost, are easy to maintain, and offer high endurance over long distances.^1–3 Various sensors can be easily integrated into underwater gliders to measure different ocean indicators, such as salinity, depth, chlorophyll, and dissolved oxygen. In recent years, with the continuous development of comprehensive ocean exploration, underwater gliders have been widely used in marine science research.⁴ However, being driven by buoyancy makes them susceptible to ocean currents, which leads gliders to often not follow arbitrary trajectories to a designated position.⁵ To reduce the impact of ocean currents on the navigation of an underwater glider, it is necessary to know the ocean current information of the glider deployment zone in advance. However, there is currently no perfect domestic ocean dynamics model for current forecasting that can be provided to gliders. Moreover, due to the power consumption and limited carrying capacity, underwater gliders are not equipped with current meters or acoustic Doppler current profilers in most cases.^6,7 To solve the above contradiction, Merckelbach et al. proposed the concept of deep averaged current velocities (DACVs), which are defined as the difference between the dead-reckoned (DR) and measured [GPS (global position system)] positions of resurfacing divided by the subsurface time.⁸ DACVs are also called deep-averaged current,⁹ deep-averaged horizontal current velocity,⁸ or deep-averaged velocities,⁵ which are the integration of ocean currents at different depths, times, and positions in an underwater glider diving profile. They can be obtained by the glider itself without relying on either digital ocean current forecasts or ocean current measuring equipment. As a unique feature of the underwater glider, the DACV forms an important component of a gliders’ scientific value. The research on DACVs can improve the operation ability of underwater gliders and promote the development of other related technologies of ocean robots.

Navigation of a glider depends on the DACVs, either to travel between waypoints or to cross currents. What’s more, the DACVs have many other applications, such as the assistance of the underwater glider path planning, current equipment calibration, and ocean model calibration. However, an accurate measurement of the DACV of a glider in real time is difficult because the accuracy depends on many factors, such as the measurement heading, pitch, repeatable hydrodynamic shape, internal waves, and biofouling.⁵ 

In many cases, the predicted DACVs are more meaningful to us than the real time DACVs. For example, only the accurately predicted DACVs can be a good aid for underwater glider navigation.^8,10 However, it is obvious that the prediction of DACVs is even more difficult. As far as the authors know, only a few scholars have tried to study the prediction of DACVs for underwater gliders and that the prediction process is complicated and can only be applied in specific sea areas. For example, Chang et al. divided the DACV into a tidal component and a nontidal component, and then, the tidal component was predicted by a tidal model, while the nontidal component was predicted by an empirical hybrid ocean model based on the experimental data collected by a glider, and the predicted DACVs consisted of the two predicted components.¹⁰ Similarly, Merckelbach divided the DACV near the northern sea into a tidal component and a nontidal component, and then, the two components were predicted by the Kalman filter based on the linear shallow water equation and the first-order Butterworth low-pass filter.¹¹ The above methods of forecasting the DACV require the assistance of sophisticated tidal models, and they only adapt to the coastal waters (where there are a lot of tidal currents). Therefore, the above reasons result in the poor universality of those DACV prediction methods. Zhou et al. regarded the DACVs as a time series and then adopted empirical mode decomposition (EMD)-based schemes to forecast the series,⁶ whereas the desired future data are assumed to be known in the decomposition process, which violates the purpose of the DACV prediction.¹² Overall, achieving the universal and effective forecasts of the DACVs of underwater gliders is a problem that cannot be solved.

In order to improve the universality and practicability of the DACV prediction method, this paper proposes a prediction method based on the analysis of the features of the DACV series, and then, the classical time series prediction method is used to forecast the DACVs. The forecast process does not require the inclusion of specific complex sea tide models, and there is no limit to the applicable sea area. Four commonly used evaluation criteria are utilized in these comparisons, and the performance of the proposed model is verified. The innovation of this paper can be concluded as follows: The features of the DACV series are innovatively identified. Furthermore, based on the features, a set of targeted methods for predicting DACVs are proposed. The research results of this paper can provide critical support for the efficient planning and accurate navigation of the underwater gliders and also provide a theoretical basis for the research on the prediction methods of ocean currents, sea winds, waves, and other ocean environmental elements matching with other types of ocean mobile observation platforms.

The remainder of this paper is organized as follows: In Sec. II, the DACVs are described and their features are indicated. A prediction framework for these features is given in Sec. III in which the prediction model and evaluation criteria are shown. The forecasting results and their analysis are shown in Sec. IV. Finally, we conclude in Sec. V.

Because underwater gliders are deployed in a zigzag motion in water, if one starts to dive from a certain position P_n, then after a diving profile, the glider will resurface at position P_n+1, which can be obtained using the GPS. Assuming that the horizontal speed of the underwater glider is known, the estimated resurfacing point P_n+1,0 in calm water can be obtained by the DR method.⁶ Suppose that T_n is the time cost of the nth profile. Then, we can calculate the DACV of the nth profile as follows:

This paper adopts four groups of DACVs derived from different gliders in recent sea trials for testing and analysis. These underwater gliders belong to a different series of sea-wing gliders manufactured by the Shenyang Institute of Automation, Chinese Academy of Sciences, the names of which correspond to SW001, SW002, SW003, and SW004. The number of DACVs of each glider is different, where SW001 corresponds to 380, SW002 corresponds to 395, SW003 corresponds to 465, and SW004 corresponds to 488. The period of the sea trial is from mid-July to the end of September 2019. The location of the sea trials is the South China Sea, which is shown in Fig. 1.

The DACV includes both the magnitude and direction; this paper uses only the magnitude for prediction. The statistical information of the DACVs is shown in Table I.

DAC is the synthesis of ocean current at different times and space positions in the profile of an underwater glider. The causes of ocean current formation are complicated,^13–15 including wind, air pressure, temperature, salinity, and tides.¹⁵ For previous research on DAC, scholars used to model DAC into three parts, namely, tidal component, non-tidal component, and wind driven component. All these components have both temporal and spatial variations, and they are rather sophisticated. The tidal component is composed of as many as hundreds of constituents. The non-tidal flows include ocean currents caused by temperature, salinity, and pressure, but the relationship remains challenging to pin down precisely.^16,17 The wind-driven component is related to sea wind, but the information of sea wind is difficult to obtain, and the wind-flow transfer function is not trivial.^18,19 It can be seen that each component of the DAC has complex causes and ingredients, so the DAC itself, of course, has more complex ingredients.

In addition to the complex ingredients, due to the limited endurance of underwater gliders, the number of successive gliding profiles is limited, so the DACV samples that are obtained by the underwater glider during a continuous voyage are also limited. For the underwater glider, most of its energy consumption occurs at the apex of zigzag motion. In order to give full play to its advantages of long-range, the depth of diving should not be too shallow in general. Nevertheless, at the same time, the deeper maximum depth means the specialized design costs. The 1000 m depth profile is the majority in specific observation tasks, and the 300 and 800 m profiles are also standard. In the 1000 m depth profile (this depth is the most common in various observation missions), the maximum number of DACVs that can be obtained will not exceed 500 in each glider deployment. Each DACV corresponds to a profile, and it takes 2–4 h to complete an observation profile (different settings, such as the pitch angle and input net buoyancy, can lead to different time costs) in most cases. Therefore, the continuous deployment time of underwater gliders will generally not exceed three months. For depths at 800, 300 m, or other, the maximum sample size of the DACVs cannot exceed 500 as well.

Although feature 1 indicates that it is difficult to quantify the pattern of DACV data, as a kind of averaged data, DACVs have apparent characteristics. In some early studies,^6,20 the DACVs were regarded as a time series, which still holds in this paper. For the time series, besides their seasonal, trend, periodicity, and other conventional characteristics, the most important is stationarity. Hence, we observe the stationarity of the four groups of DACVs, which can be judged by the augmented Dickey–Fuller (ADF) test. The ADF test is a unit root test, which can be used to judge whether a given data series is stationary or not. The existence of unit roots on a given dataset indicates an unpredictable system pattern. A more negative ADF test result means stronger rejection of the existence of a unit root for a given time series. Therefore, a negative ADF result means that the given dataset is stationary. The current generally accepted threshold is 0.05, which is also used in this study. If the length of the series is too short, the features of the series are not easy to analyze. If the length of the original series is L_org, first take the shortest length as 31, and then, the stationarity of the series consisting of the 1:L1th (L1 = 31, 32, and L_org) DACVs is determined, and the proportion of stationarity can be calculated. In the same way, we take the shortest sequence length as 32, and the stationarity of the series consisting of the 1:L2th (L2 = 32, 33, and L_org) DACVs is tested. The results of the four sets are shown in Fig. 2.

From Fig. 2, it can be seen that with the lengthening of the time series, the DACV series tends to be stable. The cut-off points of the four datasets are 189, 118, 208, and 141, respectively. After these cut-off points, all series are stable. When the above sequence shows stationarity, we further use the Ljung–Box Q-test²¹ to judge whether it is a white noise sequence, and the results show that all the above sequences are not white noises. Hence, the time series is thought to be predictable.

Random forest is an excellent machine learning algorithm, which can be used for classification and regression. In time series prediction, it has been successfully applied to many practical problems, representative applications of which can be found in many scientific fields, including engineering,²² environmental sciences,²³ financial research,²⁴ and medicine.²⁵ In these applications, small datasets are utilized, and random forest has demonstrated its excellent predictive performance. As described in Refs. 26 and 27, random forests can handle data with small sample sizes and complex component structures.

Random forest is implemented based on Bagging and decision trees. Bagging is an ensemble learning algorithm,²⁸ which randomly extracts training samples from the original sample set with replacement and trains a single weak learner. In the random forest regression model, the weak learner is a regression tree. This process is repeated to generate multiple regression trees to form a random forest, and the final prediction result is determined by the average value of the predicted values of all trees. Common decision tree algorithms include ID3, CART, and C4.5, all of which adopt the top-down greedy algorithm to construct the tree structure, select an optimal feature for splitting at each node, and recursively build the tree until the termination condition is met. In this paper, the CART algorithm is used to generate a regression tree, and the basis for selecting features and dividing points of the regression tree is the minimum mean square deviation.

The steps to build the random forest regression model are as follows:

1. Using the Bagging method to select n samples from N initial training samples (n < N) with replacement to generate m training subsets.

2. Training the regression tree using the training subset and selecting some sample features from all the sample features on the node. The left and right subtrees of the regression tree were divided according to the minimum mean square deviation, and the tree was recursively built until the termination condition is satisfied.

3. Repeat the above steps to form a random forest with multiple regression trees.

4. The samples to be predicted are input into the random forest regression model, and the average of the predicted values of all trees was taken as the final prediction result. The above steps are represented in Fig. 3.

Autoregressive integrated moving average (ARIMA) is a linearization method, which assumes the future value of a predicted variable as a linear function of past observations. Therefore, time series data that are fed to ARIMA are expected to be linear and stationary.

ARIMA performs a stationarity check on given time series data to check whether the average and autocorrelation patterns are constant over time. If the stationarity is not satisfied, ARIMA applies the difference method d times until the nonstationarity is dealt with. As a consequence, the integration order of the ARIMA model is set to d. Thereafter, the autoregressive moving average (ARMA) is applied to the result data as follows:

For the kth profile, let the actual DACV be y_k, and the random error is e_k. y_k is considered a linear function of the past p observations y_k−1, y_k−2, and y_k−p and q random errors e_k, e_k−1, and e_k−q. The corresponding ARMA equation is as follows:

In Eq. (2), the coefficient from α₁ to α_p is the autoregressive coefficient and θ₁ to θ_q is the moving average coefficient. Note that the mean of random errors e_k is zero, and the variance is constant. Similar to the parameter d, p and q coefficients are referred to the order of the model.

For the DACV series used in the prediction process, if the training sequence does not reach the critical length, it indicates that the sequence is not stationary, and then, the ARIMA model is built. If the training sequence reaches the critical length, the ARMA model is adopted directly. The determination of p and q is carried out by using the combination of the Akaike information criterion (AIC) and Bayesian information criterion (BIC) as follows:

All DACV datasets satisfy feature 1 and feature 2, in which feature 1 is targeted by RF and feature 2 is targeted by ARIMA. Considering that feature 1 and feature 2 may have certain coupling, the prediction method will also use the combination of these two methods in addition to RF and ARIMA alone. The cut-off point in Sec. II B can be used as a reasonable demarcation point, with ARIMA (RF) used for the data in front of it and RF (ARIMA) used for the data behind it. In this way, for each dataset, there are a total of four methods to predict it.

Three commonly used evaluation indicators are used to analyze the prediction results,

where x_t is the raw DACV data, $xt̂$ is the forecasted DACV data, and n is the number of DACV samples of the x_t series. Because this article implements a one-step rolling forecast, all n take 1.

We adopt the four methods proposed in Sec. III to predict the magnitude of the DACVs. Since in the practical application of underwater gliders, the persistence method (regarding the DACV in the last profile as the forecasting DACV of the next profile) is also an effective method, we choose it as a comparison model. The results are shown in Tables II–V and Figs. 4–7. It should be noted that in Figs. 4–7, in order to make the figures more clear, both the real DACVs and the predicted DACVs are sampled and drawn at an interval of 5.

From Tables II–V and from Figs. 4–7, we find that the ARIMA, RF, R-A, and A-R have higher prediction accuracy than the persistence method for all four groups of DACVs derived from different underwater gliders. The results show that the methods we proposed are effective. We can also see that ARIMA is the best for DACVs from SW001, RF is the best for DACVs from SW002, A-R is the best for DACVs from SW003, and R-A is the best for DACVs from SW004. For DACVs derived from SW001 and SW002, the hybrid model does not show an advantage, and the single model performed only slightly better than the hybrid model, which can be illustrated by the comparison between A-R and ARIMA in Table II and A-R and RF in Table III. For DACVs derived from SW003 and SW004, the hybrid model is better than the single model, which can be seen in Tables IV and V. The reason for this is that RF is well suited for predicting DACV data with feature 1, while ARIMA is highly targeted for feature 2, but there is a certain coupling between feature 1 and feature 2, which cannot be clarified. We use the position of the stationary ratio of 1 as the cut-off point, and the cut-off point of the four datasets is 189, 118, 208, and 141, respectively, which correspond to the black dotted lines in the figures. In fact, after our test, if we set the point where the stationary ratio is 0.8 or 0.9 as the cut-off point, the results remain the same. In general, the proposed method based on data statistics and the method based on machine learning are effective for DACVS data with feature 1 and feature 2.

In this paper, the features of DACVs are analyzed, and methods based on data statistics and the method based on machine learning are proposed as the targeted forecasting model. The 1000 m depth data of four different underwater gliders from recent sea trials are utilized to demonstrate the excellent forecasting performance of the proposed model. The forecasting results show that the proposed model is effective. It should be pointed out that although this paper proposes a prediction method based on the features of DACVs, these features are relatively fixed, while the methods are not fixed.

In this paper, the DACV data of the underwater glider are regarded as a time series to predict them, which can meet the vast majority of requirements of the underwater glider in scientific investigation missions. Nevertheless, the DACV data of the underwater glider are a kind of space–time data, which are influenced by human intention. For example, the operator can give the heading and set the depth in advance for the next period, which will affect the deep mean flow for the next period. Similar space–time sequence prediction can be found in Refs. 29–32. How to consider the DACV prediction of underwater gliders in this situation in more detail and realistically is the work of the next stage.

The significance of this paper is that it provides a more accurate, robust, and practical model for forecasting short-term DACVs of gliders, whereas high-accuracy DACV forecasting can be performed to improve glider navigation and path planning. The proposed model can also be applied in many other fields, such as wind speed forecasting, blood pressure forecasting, and stock price forecasting.

This work was supported by the National Natural Science Foundation of China (Grant Nos. U1709202 and 51809127), the Natural Science Foundation of Shanxi Province, China (Grant No. 201901D211248), and the Natural Science Foundation of State Key Laboratory of Robotics (Grant No. 2019-O10).

The data that support the findings of this study are available from the corresponding author upon reasonable request.