Sea state from ocean video with singular spectrum analysis and extended Kalman filter

A method for estimating key parameters of ocean waves (the dominant frequency and the significant wave height) from uncalibrated monoscopic video is proposed, based on temporal variation of the wave field, specifically time series of pixel intensities. The methodology tracks the principal component of the movement of water in the video, which we propose is associated with the dominant frequency of the ocean. To accomplish this, the singular spectrum analysis algorithm and the extended Kalman filter are used. Then, the shape of an empirical spectrum is used in order to translate the dominant frequency output into a significant wave height estimation.

This work investigates the estimation of the sea state from a single uncalibrated camera. We do not utilise techniques that are used for the prediction of the wave elevation from tidal gauges (such as harmonic analysis and wavelet networks) because the pixel intensity from video does not correspond directly to wave elevation. The pixel intensity can be considered proportional to the lights reflected from the water surface [11].
Remote sensing from simple video cameras has been widely applied for acquiring the nearshore hydrodynamics and morphology [7,14,26]. The bathymetry is estimated with a celerity-based depth-inversion method that utilises the dispersion relation of shallow water and the spatial correlation of pixel intensity signals indicating propagation of waves. Based on this information, the nearshore sea levels [9,18] and current predictions [22] are acquired.
The present work presents a technique that is applicable to real environments (unlike [11]), deep water (unlike [18,19,27,29,30]), does not use in situ devices for calibration (unlike [11,30]) and is validated with videos that have corresponding in situ measurements in a variety of sea states (unlike [19,[23][24][25]). Unlike the techniques that estimate hydrodynamics, morphology or sea state in the nearshore, this work does not utilise information of foam from breaking waves. (In deep water, foam is present in very high sea states.) In our previous work [16], we use the linear Kalman filter and the least squares approximate solution in order to form the uncalibrated ocean video amplitude spectrum. We then use ocean theory in order to calibrate this spectrum into metres and estimate the significant wave height. In the present work, the goal is to solve the same problem with a novel methodology. We verify the sea state estimations with the same video data as before.
We start, in Sect. 2, by providing an introduction to key ocean theory used by the methodology, before describing the methodology itself (Sect. 3), and then demonstrating its efficacy on real videos of the sea, with in situ buoy measurements for validation (Sect. 4).

Ocean theory: the Pierson-Moskowitz spectrum
The Pierson-Moskowitz spectrum [21] is an empirical spectrum of the ocean formed from data acquired from accelerometers installed on weather ships. The spectral energy in terms of angular frequency ω is expressed as: where α = 8.1 × 10 −3 , β = 0.74, g is the gravitation acceleration and ω 0 = g/U , where U is the wind speed at 19.5 m above the ocean surface. The dominant angular frequency ω m is equal to: The area under the spectrum is equal to the integral of the function: The significant wave height can be found to be equal to four times the square root of the area under the spectral density [3].

Methodology
The aim of this work is to track the main oscillatory component from video time series of pixel intensities that is associated with the ocean's movement. This enables the estimation of the ocean dominant frequency and the significant wave height. To achieve this, a methodology is introduced that combines the SSA algorithm and the nonlinear Kalman filter. It also incorporates ocean theory presented in Sect. 2. Fig. 1 Example of matrix of w-correlations of SSA algorithm from ocean video. From this example, a correlation of certain PCs is visible, such as PCs 1-2, 3-4, 5-6, 7-8. Additionally, some higher PCs can be considered to contain more noise (PCs 11-max)

Singular spectrum analysis (SSA) algorithm
Historically, the SSA algorithm is associated with work published in the 1980s, e.g. [10]. In the context of time series analysis, the SSA algorithm decomposes the input signal into a set of additive components, which are labelled as either trend, oscillatory or noise components.
In the context of this work, time series of pixel intensities is given to the SSA algorithm. The first four elementary reconstructed components (RCs) are summed to provide a new time series, which is given to the extended Kalman filter in order to estimate the dominant frequency. The hypothesis is that the SSA algorithm will concentrate the information of the central component from the video, which is associated with the dominant wave of the ocean, in the first RCs.
Although it would be expected that the dominant frequency is isolated in the first RC, practically this was not found to be the case. Empirically, selecting the first four in all cases (see Sect. 4.2) was sufficient. This selection is also based on observations from the matrix of w-correlations (see Fig. 1). Specifically, in many cases it was observed that the first four RCs had a strong correlation.
The next step involves the determination of the dominant frequency from the sum of RCs 1-4 of the SSA algorithm. Practically, the Fourier transform of this time series includes more than one peak. Selecting the highest peak does not in all cases correspond to the dominant frequency of the ocean (see Fig. 2). Since this selection of RCs usually includes more than one wave, determining one frequency is not a straightforward task. This is the reason the extended Kalman filter is used in the following step.

Extended Kalman filter algorithm
The extended Kalman filter algorithm with the environment definition described in the following text is very efficient at distinguishing one main frequency and isolating the remaining video elements as noise. As mentioned in the previous section, the goal is to identify the principal component of water movement from the sum of RCs 1-4, which is hypothesised to be related to the dominant wave of the ocean. The true signal is a sinusoid: where a is the amplitude, φ the phase, and t is the time. The derivative of the signal with respect to time is equal to: and the second derivative: The second derivative can be expressed as a function of the angular frequency and the true signal as: This is useful in the context of this work, because it does not include amplitude and phase, and instead focuses on the true signal. Additionally, the sinusoidal form of the signal is now included in the environment definition. The environment definition is: The Jacobian of matrix g is computed at each time step with the current estimates of the states. By taking the partial derivatives of the matrix, the system's dynamics matrix F is found to be equal to: whereω is the predicted value or estimate of ω and similarlŷ x is the predicted value or estimate of x at the current time step. The fundamental matrix is not used for propagating the states, but rather only for the calculation of the gains, and is approximated with the first two terms of Taylor series. With the environment definition described here, the nonlinear Kalman filter can be solved as in [2]. The derivative of the angular frequency is equal to zero. The true state of the angular frequency is constant because the sea state is not expected to change in the duration of the video. Although the unknown true state of the angular frequency is constant, the algorithm's estimation of this value varies at each iteration, as can be seen in Fig. 3c. Including the angular frequency is our environment definition is important because the value estimate is used for inferring the value of the significant wave height in the next step of the methodology.
The Kalman filter outputs a value of the unknown angular frequency. This angular frequency can be used directly as the ocean dominant angular frequency. In the following section, this value is given as input to the theory of the Pierson-Moskowitz spectrum in order to get a value of the significant wave height. As a side note, the described methodology is performed on one pixel time series. For acquiring more accurate and reliable dominant frequency estimations, a set of pixels (or all pixels) can be used individually and the dominant frequency is found as the average. In the case of the experimental results of this work (see Sect. 4), a set of pixels equal to the image width is used. Figure 3 demonstrates the functionality of methodology. The time series of pixel intensities is used as input to the SSA algorithm, and the principal component of the movement of water in the video is speculated to be included in the sum of RCs 1-4 (Fig. 3a). This main component is isolated from all other video components with the extended Kalman filter (Fig. 3b). The filter provides an estimate of what we hypothesise to be the ocean dominant angular frequency in each time step (Fig. 3c) and the limits of certainty for that prediction from the square root of the corresponding diagonal element of the covariance matrix (Fig. 3d).

Significant wave height
The significant wave height (h s ) can be found as equal to four times the square root of the area under the ocean spectral density [3]. The dominant frequency from the previous steps can be given as input to the Pierson-Moskowitz spectrum equation (see Sect. 2). The shape of this empirical spectrum can then be used in order to approximate h s . From Eqs. (2) and (3), h s can be expressed in terms of dominant angular frequency ω m (found in the previous steps) as: where α = 8.1×10 −3 and β = 0.74. The values of dominant frequency and significant wave height are the outputs of the methodology.

Experimental results
Two sets of video data that have corresponding in situ measurements are used in order to test the accuracy of the proposed technique. Both sets comprise videos with duration of approximately one minute. The first set is taken from a shipborne camera in experiments done on the 24 November 2014 in the North Atlantic sea. Two buoys measured the ocean at the same time in a nearby location. The sea state in this set of videos is approximately the same, as the state is not expected to change in a large degree in the time span of a few hours. The significant wave height of the shipborne video is approximately 3.1m-3.4m. The second set of video data is taken from a camera on the Frying Pan Shoals tower, a former lighthouse located approximately 39 miles southeast of Southport, North Carolina. A 24-h live video footage of the ocean is available online [13]. Although the camera is panning showing a panoramic view,

Preprocessing
The preprocessing step involves stabilisation. For the shipborne video, the rotational movement of the ship (pitch, roll, yaw) is stabilised by stabilising the horizon. This is achieved with the rotational tracker of Adobe After Effects [1]. Two rectangles are drawn above the video in order to stop the tracking points from moving horizontally while tracking the horizon. Figure 4a shows a typical frame from ship video and the rectangles drawn in order to track the horizon, and Fig. 4b presents how each tracking point is selected. For the tower video, the video is stabilised only in cases of high local wind with the stabilisation features of Adobe After Effects. A single set of pixel locations is used for computational efficiency as in Fig. 4c and f. Figure 4c presents the result from ship video after preprocessing. Figure 4d shows a typical frame from tower video and Fig. 4f the result after preprocessing.

Main methodology
The first set of data from the ship examine the behaviour of the methodology for an approximately statistically stationary sea state. The second set of data from the tower examine the behaviour for a variety of sea states, as the videos were captured in different days. The SSA algorithm is run in all cases for a window length of 350, which is determined empirically.
The matrix of w-correlations provides a good indication of whether the window is too small or too large. Specifically, in cases of smaller window size the association between RCs in the main diagonal is weak (the model is too general). With a larger window size, high values are concentrated in positions further away from the main diagonal. In this case, it can be interpreted as the algorithm is overfitting. From empirical observations, the estimations are relatively insensitive to the window length. That is, even if different window lengths are used (for example 250 or 450) the impact on the estimations is minimal (see Sect. 4.3).
The shipborne results are presented in Fig. 5. Each point represents a one minute video captured in the same day. The buoys were deployed on the sea surface after 9:15 a.m., any videos before then are presented here only to show the behaviour of the methodology. As mentioned, the sea state is not expected to change in a large degree in the span of 1 h. The error metrics between the video estimations and buoy measurements are presented in Table 1.
From Fig. 5, it is observed that the video estimations vary in values but are close in proximity to what the buoys indicate to be the true sea state. From Table 1, the 0.19 m of RMSE and 4.83% of MAPE indicate that the video estimations are not very distant from the buoy measurements. Until this point, the results support the hypothesis of the present work that   the methodology estimates the significant wave height in a satisfactory degree of accuracy. The tower video are used for examining the behaviour of the method across a variety of sea states. The tower video results are presented in Fig. 6. Each point represents a video captured on a specific date. The technique estimates lower dominant frequencies (higher h s ) for higher sea states and higher dominant frequencies (lower h s ) for lower sea states, as expected.
The error metrics from tower video are presented in Table 2. Although the error metrics of 0.23 m RMSE and 15.29% of MAPE are higher than the ones observed from shipborne video, they still remain in acceptable levels in showing that the estimations are meaningful. It should be mentioned that the MAPE metric provides higher error when the pairs of true-estimated values are lower. This is one possible reason for the higher value of MAPE with the tower video, as the tower video is captured in both lower and higher sea states. The experimental results show the methodology's potential for estimating the significant wave height.

Sensitivity analysis
From the matrix of w-correlations (see Fig. 1), it is observed that in many cases the first two or the first four elementary reconstructed components (RCs) have a strong correlation. It is also observed that components after RC 10 contain more noise. Additionally, there is a strong correlation between different components, such as RCs 3-4 and RCs 5-6. From empirical observations with multiple videos, the selection of only the first two RCs does not provide accurate estimations. The selection of RCs 1-4, 1-6 and 1-8 provides more accurate estimations. From RCs 1-10 and higher the estimations become less accurate. Possible reason for this effect is the inclusion of more noise as the component number increases. For example, see Fig. 7.
In terms of the window length parameter of the SSA algorithm, values have been tested ranging from 50 to 1000 (the number of frames is approximately 1000 for the shipborne video and 2000 for the tower video). For values between 100 and 850, the sea state estimation does not vary in a large degree. For example, see Fig. 8.
The speculation for the less accurate results with a small window size is that the model is too general; that is, the association between the different elementary reconstructed components (RCs) is not clearly defined. Similarly, the speculation for the less accurate results with a very high window size is that the model is too specific; that is, the association between the different RCs is overspecified. Fig. 7 Example of sensitivity analysis with a varying selection of RCs of the SSA algorithm with tower videos as input. The root mean square error (RMSE) for the estimation of the significant wave height Fig. 8 Example of sensitivity analysis with a varying window size of the SSA algorithm with tower videos as input. The root mean square error (RMSE) for the estimation of the significant wave height A two-term Taylor series is used for approximating the fundamental matrix, which is only used for the calculation of the Kalman gains. The sea state estimation is approximately the same with the use of higher-order Taylor series. Up to five-term Taylor series have been included for the sensitivity analysis, with no significant improvements in the estimation of the sea state. Including more terms could potentially be beneficial if the approximate fundamental matrix was used for propagating the states. However, in this work the propagation of the states is achieved with integration of the differential equations over the sampling interval.

Conclusion
The range of significant wave height for which the technique is tested is between 0.5 and 3.6 m. Higher sea states are not examined, and thus, the applicability in such instances is not known. Testing the technique in stable and varying sea states shows the potential of the methodology at estimating the sea state from video.
Future work will include further testing for higher sea states. In the context of practical utilisation of the present work, in some cases estimation of higher sea states might not be required. For example, in the case of use of the work for the execution of maritime operations, very high sea states can already be identified.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.