Extracting clearer tsunami currents from shipborne Automatic Identification System data using ship yaw and equation of ship response

We have explored tsunami current signals in maritime Automatic Identification System (AIS) data during the 2011 Tohoku, Japan, tsunami. The AIS data were investigated in detail taking into account ship motion and response to tsunami current. Ship velocity derived from AIS data was divided into two components in terms of the ship heading: heading-normal and heading-parallel directions. The heading-normal velocity showed good agreement with the simulated tsunami current, as mentioned in our former research. Here, we found the heading-normal velocity was contaminated by non-tsunami noises that were mostly related to the ship yaw motion around the pivot point. The noises due to the yaw motion were reasonably corrected in the heading-normal velocity. The corrected heading-normal velocity clearly showed better agreement with the simulated tsunami current. Although the heading-parallel velocity is basically the navigation speed, and is mostly controlled by ships’ captain, we could find the heading-parallel velocity was also drifted by tsunami currents. The corrected heading-normal velocity was still a ship response to the tsunami current. Based on an equation of a ship response to tsunami currents, we numerically estimated tsunami current from the corrected heading-normal velocity. We could find very slight improvements in estimating the tsunami currents, which indicated that this operation possibly worked as a secondary correction. Tsunami currents of tens of centimeters per second are expected to be suitably detected using AIS based on discussion on detection limit.


Introduction
There are many kinds of observations of tsunamis. We have mostly measured sea surface elevation for quantitative evaluations (Satake 2007;Joseph 2011). Measuring sea surface height is fundamental to understand and utilize the oceans. Many practical methods have been developed and utilized to measure sea levels including tsunamis (IOC 2006). Measuring water surface levels in coastal tide wells has a long history (Satake et al. 1988;Rabinovich et al. 2006;Namegaya et al. 2009). Aerial acoustic waves and/or frequency-modulated continuous waves (FMCW) reflected at the sea surface are used to detect water surface location (Míguez et al. 2005;Joseph 2011;Woodworth 2016). A floating buoy with a Global Navigation Satellite System (GNSS) receiver has been able to measure water surface position (Chadwell and Bock 2001;Kawai et al. 2013;Kato et al. 2018). A coastal GNSS station can detect water surface with coherent direct and reflected rays (Anderson 2000; Löfgren and Haas 2014;Larson et al. 2017). Methods from below sea level and sea bottom are utilized as well. Bubble gauges were used to estimate water column height based on hydrostatic pressure equilibrium (Pugh 1972;Woodworth and Smith 2003). Travel time of ultrasonic wave Inazu et al. Earth, Planets and Space (2020) 72:41 reflection is used to measure the distance between sea bottom and sea surface (Takayama et al. 1994;Hashimoto et al. 1996). Pressure gauges at the deep-sea bottom feel the total water column height above the seafloor (Filloux 1983;Eble and Gonzalez 1991;Hino et al. 2013). There are satellite communication systems such as the Deepocean Assessment and Reporting of Tsunamis (DART) system (Bernard et al. 2006;Mungov et al. 2013), and seafloor cable systems such as NEPTUNE Canada, DONET, and S-net (Rabinovich and Eblé 2015;Kaneda et al. 2015;Kanazawa et al. 2016). These are based on fixed-point observations.
Meanwhile, moving-point observations have been increasingly considered as well. Satellite altimetry has actually detected tsunamis (Okal et al. 1999;Godin et al. 2009;Song et al. 2012). Aircraft radar altimetry will be able to detect tsunamis (Wright et al. 2009;Hirobe et al. 2019;Mulia et al. 2020). GNSS-reflectometry using flying vehicles will be used to measure sea surface heights including tsunamis. Reflection receivers widely include low Earth orbit (LEO) satellites (Clarizia et al. 2016;Wickert et al. 2016;Mashburn et al. 2018), manned aircraft (Lowe et al. 2002;Ruffini et al. 2004;Cardellach et al. 2014), and unmanned aerial vehicles (Ichikawa et al. 2019). GNSS attached to a seagoing ship could detect a 10-cm tsunami height (Foster et al. 2012) because the GNSS positioning involved ~ 10-cm precision (Foster et al. 2009). If ordinary GNSS facilities with meter-order precision attached to present commercial ships were replaced by such high-precision GNSS facilities, geophysical signals as well as tsunamis would be detectable from many unspecified ships (Roggenbuck et al. 2014;Inazu et al. 2016).
Tsunamis involve water currents mostly in the horizontal direction. The water currents due to tsunamis have been measured by acoustic Doppler current profilers (ADCPs) at fixed points (Bricker et al. 2007;Lacy et al. 2012;Sobarzo et al. 2012). Terrestrial high-frequency (HF) radar has been used to spatially cover horizontal water currents within tens of kilometers from the radar stations (Hinata et al. 2011;Benjamin et al. 2016;Dzvonkovskaya et al. 2018). These facilities are mostly deployed for the purpose of monitoring coastal ocean circulations with wide applications (Roarty et al. 2019). Inazu et al. (2018) found that tsunami currents could be detected using navigation data derived from the Automatic Identification System (AIS). They demonstrated a feasibility of tsunami detectors from crowd-sourced moving vessels. If shipborne ADCP data are widely shared, ocean currents including tsunamis will be suitably detected at moving vessels (Horiuchi et al. 2015;Miyazawa et al. 2015). However, according to the current maritime convention (IMO 2016), AIS does not include ADCP facilities and data. Such data can be only aggregated from specified vessels. On the other hand, AIS data are aggregated from a large number of unspecified vessels thanks to regulations of the International Maritime Organization (IMO 2016). We have investigated this versatile AIS for tsunami current measurements (Inazu et al. 2018). In the present study, a more efficient method incorporating ship motions is demonstrated to detect clearer tsunami currents.

AIS data
We use AIS data during the 2011 Tohoku, Japan, tsunami, which are the same data as used in the authors' former research (Inazu et al. 2018). Based on the International Convention for the Safety of Life at Sea (SOLAS), vessels whose gross tonnages are basically greater than 300 gross tonnages have obligation to transmit AIS messages using very-high-frequency (VHF) maritime bands (IMO 2002(IMO , 2016. Transmitted AIS messages can be received by terrestrial stations with < ~ 100 km distances, and be received by LEO satellites (Høye et al. 2008;Eriksen et al. 2010;Carson-Jackson 2012). In March 2011, around Japan, there were only terrestrial AIS data available. Temporal sampling of transmitted AIS messages is typically tens of seconds during navigation. A formatted sentence of AIS includes dynamic and static information. Dynamic information includes time, horizontal position (latitude/longitude), speed over ground (SOG), course over ground (COG), ship heading (HDG), and so on. Static information includes Maritime Mobile Service Identity (MMSI), ship type, ship name, and so on. The requirement for the AIS format was given by the IMO Maritime Safety Committee (IMO 1998). Specifications are defined in documents of the International Telecommunication Union (ITU) and the International Electrotechnical Commission (IEC). They have been sometimes updated, and recent versions are shown in ITU (2014) and IEC (2018). Brief descriptions are also available (Tetley and Calcutt 2001;Masuda et al. 2014;IMO 2016). There are almost real-time ship distributions derived from AIS data available on websites such as MarineTraffic (https ://www. marin etraffi c.com), VesselFinder (https ://www.vesse lfind er.com), Shipfinder (http://jp.shipfi nder .com/Monit or), and so on.
Regarding the precision of position of the AIS transponder, Differential GNSS (DGNSS) correction (i.e., Msg. 17) is suitably applied when the AIS transponder can find an available reference station. Ordinary GNSS single-point positioning is applied if there is no available station (IEC 2018). When DGNSS correction is applied, positioning precision is expected to be < 10 m. Resolution of the heading direction derived from compass seems not to be strictly determined (IEC 2018). In most cases in AIS data, HDG and COG are obtained in 1.0 and 0.1 degrees, respectively. Figure 1 shows the ship distribution during the 2011 Tohoku tsunami derived from the AIS data. All the identified vessels were > 300 gross tonnages. Smaller vessels such as fishery boats were not included in the AIS data although there were probably a great number of such smaller vessels. At offshore regions of Iwate and Fukushima, AIS data could be retrieved at terrestrial stations during < 45 min after the earthquake generation. Meanwhile, offshore of Miyagi, AIS data received at terrestrial stations were lost immediately (a few minutes) after the earthquake generation probably because of power failure due to the strongest ground motion near the epicenter. As shown in Inazu et al. (2018), we found 16 ships navigating offshore Tohoku (depth > ~ 100 m), and clear tsunami current signals from 13 ships (#1-#13) ( Table 1). These ships actually experienced tsunami currents far greater than 0.1 m/s, as also indicated by tsunami simulation (Fig. 1).
The AIS data are investigated in detail (Fig. 2). SOG, COG and HDG are obtained from AIS records in addition to the time and position. In principle, SOG and COG are calculated from the time derivative of the ship positions. HDG is obtained by gyro/ magnetic compass. The rate of turn (ROT) or angular velocity of the ship is obtained from the time derivative of HDG. The HDG-normal velocity (sway v n = SOG sin(COG − HDG) ) is calculated as well as a b Fig. 1 Ship distributions 24 min after the 2011 Tohoku tsunami generation identified from AIS data, and simulated maximum tsunami current. The rectangular area in a is enlarged in b. Labeled 13 ships (#1-#13) are used to evaluate tsunami currents in this study (Table 1). Another three ships (*1-*3) are not used due to data quality issues (Inazu et al. 2018). Dashed isobaths of 100 m and 500 m are shown in addition to the Japan Trench.
Solid contours indicate the cumulative tsunami source of Hossen et al. (2015) Inazu et al. (2018). The HDG-parallel velocity (surge v p = SOG cos(COG − HDG) ) is also calculated. Figure 3 shows ship velocity and simulated tsunami current in the HDG-normal direction. Initial condition or generation process of the tsunami is a crucial factor for the tsunami simulation result. Simulation results using two initial conditions are compared. One is given by Saito et al. (2011) which estimated an instant tsunami generation. The result was used in Inazu et al. (2018). The other is given by Hossen et al. (2015) which estimated a dynamic tsunami generation during 3 min. Both simulation models well represent v n . Saito et al. (2011) fail to correlate #1 but correlate well #5 (Fig. 3a). Meanwhile Hossen et al. (2015) correlate well #1 but fail to correlate #5 (Fig. 3b). Overall correlation coefficients are equivalent as much as 0.72-0.73 derived from both the models. Both the models were previously compared by Baba et al. (2017), and the model of Hossen et al. (2015) was able to show possibly better waveforms in both near-and farfield stations in the Pacific Ocean. In the present study, the simulation result derived from Hossen et al. (2015) is mainly used to examine the AIS data. The simulation result derived from Saito et al. (2011) is complementally used to evaluate sensitivities of the analyses (Additional file 1: Figs. S2 and S3).
The AIS data are shown with the tsunami simulation in Fig. 4 and Additional file 1: Fig. S1. The AIS data shown are averaged over a 1-min time window, and resampled by 0.25 min. As shown in Fig. 3, the HDG-normal velocity (v n ) is a good proxy of tsunami currents in this direction. According to theoretical considerations by Inazu et al. (2018), it is reasonable that the HDG-normal velocity immediately (< a few minutes) responds to longwave tsunami currents with equivalent speeds. Here, the HDG-parallel velocity (v p ) is evaluated as well. The HDGparallel velocity is basically recognized as forward speed of ship (Fig. 4e, f ). According to the major tsunami warning which was first issued 3 min and next updated 28 min after the earthquake generation (Tsushima et al. 2011;Ozaki 2012), captains of ships had to carefully judge how the tsunami came. Several ships perhaps adjusted to  decrease their speeds (#3, #6, #7, #10, and #12) after the first tsunami warning. Several ships (#3, #11, and #12) also decreased their speeds around the next update of the tsunami warning ( Fig. 4 and Additional file 1: Fig. S1). Thus, captains of these ships probably controlled their speeds or the HDG-parallel velocities. It is difficult to distinguish the tsunami and man-made control in the HDGparallel velocity. Nevertheless, ships such as #1, #2, and #5, where man-made controls seem not drastic, clearly showed that the HDG-parallel velocities were drifted by tsunami currents with equivalent amplitudes (Fig. 4e and Additional file 1: Fig. S1). This indicates that the response of the HDG-parallel velocity to the tsunami currents is qualitatively similar to that of the HDG-normal velocity.

Rate of turn and heading-normal velocity
There are still non-tsunami deviations in the HDG-normal velocity (e.g., #3 and #11 in Fig. 4d). When the HDGnormal velocity is used as a tsunami current meter, such deviations will be significant noises. Here, we considered ROT derived from the time derivative of HDG, and found clear (negative) correlations between ROT and the HDG-normal velocity (Fig. 4c, d).
We try to understand the relationship of the negative correlation employing a simplified model (Fig. 5) in which a vessel moving forward steers rightward in a calm ocean. First, the vessel is moving forward. HDG should be the same as COG (Fig. 5a). Second, the vessel is abruptly steering rightward or yawing clockwise. When HDG is increasing, ROT is positive. Probably COG is left behind HDG just after steering (Fig. 5b). Thereafter COG catches up with HDG (Fig. 5c). While the pivot point of a stopped ship is located at the center of mass of the ship, the pivot point of a ship moving forward is moving forward as well (Tzeng 1998;Clark 2005;Seo and Mishu 2011). In addition, a GNSS antenna is usually attached near the navigation bridge, and the bridge is located at the back for many ships, in particular for ordinary commercial ships such as cargo ships and tankers (Chihara 1969). Ship position and COG are strictly related to the location of the GNSS antenna. Thus, the rotation axis of HDG (i.e., pivot point) and the rotation axis of COG are separated by a distance comparable to the ship length. So, there is an apparent sway motion during yawing (Fig. 5). The apparent sway speed is where L is the ship length. A is a coefficient which seems mostly to be 1.0 but may be variable. When ROT is positive (i.e., HDG is increasing), v s indicates negative direction of the heading-normal velocity (Figs. 2 and 5b). It is known that, when an advancing ship turns its heading, the pivot point moves forward to keep the circle radius of the pivot-point trajectory effectively shorter than that of the center of mass trajectory (Seo and Mishu 2011). The essence of the apparent sway has been known in turning circle tests in naval architectures (Kijima et al. ( We find good agreement between the HDG-normal velocity (v n ) and the apparent sway velocity (v s ) in particular before tsunami arrival times (Fig. 6a). The components derived from ROT are well corrected from the HDG-normal velocity, in particular for several ships of #3, #4, #6, and #11 (Fig. 6b). The corrected HDG-normal velocity (v = v n − v s ) shows better agreement with the simulated tsunami currents (v c ). Correlation coefficients between v n − v s and v c (mean value of 0.78) become greater than those between v n and v c , (0.72) for almost (11/13) of the ships (Figs. 3b and 6b). Such improvements in the correlation coefficients are also found for 10/13 ships when tsunami simulation employed the initial condition of Saito et al. (2011) (Fig. 3a and Additional file 1: Fig. S2). The mean correlation coefficient similarly increases to 0.78 from 0.73. These indicate that the correction using ROT robustly removes non-tsunami components in the AIS data.  3). The correction is expected to be relatively effective in order to find smaller tsunami signals. Thus, the HDG-normal velocity could be separated into two components: one is the ship response to external force (tsunami current), and the other is the apparent sway velocity due to ship yawing.
Here, A ~ 1.0 is roughly employed. However, as mentioned above, this coefficient is essentially variable, depending on forward speed, ROT, and so on. Although it was noted above that the GNSS antenna with the navigation bridge is located at the back for many commercial ships, the position of the GNSS antenna actually depends on respective ships. So, A should be adaptively determined as a temporally changing variable for each vessel. It may need further work to propose a method to adaptively obtain the coefficient using time series with a certain time length.

Simulating ship response and tsunami current
The corrected HDG-normal velocity is still not certainly the tsunami current in this direction (Fig. 6b). The corrected HDG-normal velocity should be exactly the ship response to the tsunami current. An equation of motion of a floating body was derived to examine an analytical response of the ship (Inazu et al. 2018). Here, the equation of motion is again used to estimate the ship response to the tsunami current in the HDG-normal direction:

Fig. 6
Ship velocity and tsunami current in the HDG-normal direction. In a, green and black lines indicate the HDG-normal velocity shown in Fig. 3 (v n ) and the apparent sway velocity due to ship yawing in Fig. 5 (v s ), respectively. In b, green and red lines indicate the corrected HDG-normal velocity (v = v n − v s ) and the simulated tsunami current (v c ) shown in Fig. 3b, respectively. Arrows indicate periods where v s is relatively large in a and is effectively corrected in b. Correlation coefficients between v and v c are attached in the same way as in Fig. 3 where v and v c are ship response velocity and tsunami current, respectively. C D and C M are constant coefficients of drag force and inertia force (Dean and Dalrymple 1985;Reddy and Swamidas 2013). B is the breadth of the ship. As described in the equations (2)-(4) in Inazu et al. (2018), the ship mass is proportional to the product of length, draft, and breadth of the ship, and the drag force in the HDG-normal direction is proportional to the product of length and draft of the ship. Thus, the ship acceleration in this direction is only dependent on the ship breadth, or proportional to 1/B. The differential Eq.
(2) is implemented using a fourthorder Runge-Kutta method to numerically calculate v and v c . When v c is given, v is simulated. Conversely, when v is given, v c can be simulated as well. Figure 7 shows a result of a test numerical experiment. Parameters are assumed as follows: C D = 1, C M = 2, and B = 20 m. First, tsunami current (v c ) is given as an arbitrary input. The ship response (v′) can be simulated. There is a reasonable ship response (v′) to a steep tsunami current (v c ) (Fig. 7a). While v c is forced by arbitrarily changes, v c − v′ reduces to be asymptotically zero (Fig. 7b). The response of v′ is understandable as indicated by Fig. 6 of Inazu et al. (2018). Next, the tsunami current (v c ′) can be conversely simulated from the simulated ship response (v′). Prime marks (e.g., v c ′ and v′) indicate simulations. v c and v c ′ should be precisely consistent according to the analytical solution. In the test experiment, v c includes both abrupt (times at 35 and 45 min) and moderate (other times) changes. Although there are numerical errors between v c and v c ′ probably due to the used differential method, v c ′ mostly agrees v c during the moderate changes in v c . The errors tend to be distinct during the abrupt changes in v c . Even the fourthorder differentiation has difficulty in perfectly representing such abrupt changes that may be suitably represented by high-order components.
This method is applied to the data in Fig. 6b. Two simulations are carried out. One is that the tsunami current (v c ) is used to simulate the ship response (v′) (Fig. 8a). The other is that the corrected HDG-normal velocity (v = v n − v s ) is used to simulate the tsunami current (v c ′) (Fig. 8b). In a practical application, we expect that v is calculated from AIS data, and v c ′ is estimated from v. B is taken from each ship (Table 1). Simply, C D = 1 and C M = 2 are assumed.
We see estimated ship velocity (v′) in the HDG-normal direction responds to tsunami current (v c ) with a response time less than a few minutes (Fig. 8a). Estimated tsunami current (v c ′) to force the ship velocity (v) seems to be reasonable as well (Fig. 8b). Then, ideal results are v c ′ = v c , and v′ = v. Correlation coefficients between v c ′ and v c are found to be better than those between v and v c for 7/13 ships, and the mean correlation coefficient may be slightly improved to 0.80 from 0.78, but these improvements are not so significant (Figs. 6 and 8b). Correlation coefficients between v′ and v were better than v and v c for 8/13 ships, and the mean correlation coefficient is as much as 0.78 (Figs. 6 and 8a). When v c is derived from Saito et al. (2011), the mean correlation coefficients between v c ′ and v c and between v′ and v both perversely get worse to 0.71 and 0.70, respectively (Additional file 1: Figs. S2 and S3). The differences in the correlation coefficients between Fig. 8 and Additional file 1: Fig. S3 are small in terms of the waveforms but are attributed to slight phase deviations (< a few minutes) which are probably caused by the difference between the initial conditions of the two models. A dynamic generation process during 3 min was carefully involved in the initial tsunami condition by Hossen et al. (2015). Only the tsunami simulation using the dynamic initial condition possibly yielded slightly better representations of v c ′ (and v′) (Fig. 8). Perhaps, the simulation considering the ship response by Eq. (2) potentially works to estimate tsunami current from AIS data. However, to clearly confirm that other conditions should be simultaneously improved.
Following factors that are uncertain and/or incomplete should be also considered to better represent v c ′ and v′. Numerical errors in the numerical simulations for the ship response may not always be negligible in particular for discontinuous changes of the tsunami (Fig. 7). a b Fig. 7 Test numerical simulation based on Eq. (2). In a, red, green, and blue lines indicate tsunami current given arbitrarily (v c ), ship response (v′) simulated from v c , and tsunami current (v c ′) simulated from v′, respectively. In b, red and blue indicate v c -v′ and v c ′-v′, respectively Inazu et al. Earth, Planets and Space (2020) 72:41 The equation of motion used was useful to understand essential behaviors between tsunami current and ship response, but the equation may be too simplified. Equation (2) may be better to take into account a total velocity that simultaneously considers both the HDG-normal and HDG-parallel velocities, because of non-linearity of the equation. These may modify C D and C M (Sarpkaya and Isaacson 1981;Dean and Dalrymple 1985;Reddy and Swamidas 2013). Equations for ship maneuvering (Yasukawa and Yoshimura 2015; Perara 2017) may be useful to rigorously simulate ship behaviors. However, they need a number of parameters, and much of them are not available in AIS data. At this moment, we note that the correction using the yawing effect is first necessary to obtain clearer tsunami current signals from AIS data (Fig. 6). Although the estimation of the tsunami current from the corrected HDG-normal velocity using Eq. (2) is also reasonable in principle, this operation yielded a slight improvement, and then will be a secondary correction (Fig. 8b). Other factors mentioned above including accuracy of tsunami source probably involve certain amounts of uncertainties in estimating the tsunami current using the equation of ship motion.

Detectability of tsunami currents from AIS data
We have clearly found tsunami currents greater than ~ 1.0 m/s in the HDG-normal direction (Figs. 4 and 6). The detectability of tsunami currents from AIS data is discussed. According to Figs. 4 and Additional file 1: Fig.  S1, when |COG − HDG| > 5°, it is quite extraordinary. SOG is typically 5-10 m/s.
To examine the detectability of the tsunami whose time scales are mostly less than 1 h, we calculated standard a b deviations (sigma) of COG-HDG during 2 h before the earthquake, and they were 0.3-1.0° (Table 1). Based on a three-sigma limit, when |COG − HDG| is greater than 1.0-3.0°, they should be fairly extraordinary. Similarly, the three-sigma limits of v = v n − v s were 0.1-0.3 m/s, which should be lower limits of tsunami current detection using current AIS data (Table 1). According to the small-amplitude water-wave theory, the relationship between horizontal water velocity and tsunami height is well-known when the long-wave approximation is assumed (Dean and Dalrymple 1985): where v is long-wave tsunami current speed, η tsunami height, h water depth, and g gravity acceleration.
When v = 0.1-0.3 m/s and h = 100 m are assumed with g = 9.8 m/s 2 , an equivalent tsunami height is η = 0.3-1.0 m. Once offshore tsunami height is estimated, coastal tsunami height is shortly calculated using a kind of Green's law. According to experiences of offshore and coastal tsunami height observations for several earthquakes (Hayashi 2010), we note that tsunami currents due to near-field megathrust earthquakes whose moment magnitudes are greater than about 8 will be suitably detected by current AIS data. A survey by Reynolds (2012) reported that a container ship in offshore navigation showed |COG − HDG| < 5° in most cases during days. Our analysis of the three-sigma limits of |COG − HDG| was 1.0-3.0°, which is roughly smaller than 5°. This is reasonable since our analysis was carried out using only 2-h time series for the tsunami detection.

Summary
The AIS data during the 2011 Tohoku tsunami were investigated in detail after Inazu et al. (2018). Ship velocity during offshore navigation strongly deviated with mostly equivalent tsunami current velocity not only for the heading-normal direction but also for the headingparallel direction. The heading-normal velocity of the ship was affected by the heading rotation because the pivot point of the heading and the GNSS antenna are separately located with a distance comparable to the ship length. This effect was suitably corrected, and the corrected heading-normal velocity showed better agreement with simulated tsunami currents. However, the corrected heading-normal velocity was still not tsunami current itself, and should be exactly the ship response to the tsunami current. Based on a theoretical equation of motion of a floating vessel, we could numerically estimate tsunami current from the corrected heading-normal velocity. Improvements were very slight or not so significant in correlation between the estimated tsunami current and tsunami wave simulation, probably because other uncertainties including tsunami source model might affected in estimating the tsunami current. At this moment, the heading-normal velocity derived from AIS data will be more suitably used for tsunami current meters once the yawing effect is corrected. We expect to find tsunami currents of tens of centimeters per second using AIS, and hope to examine smaller tsunamis.