SEA STATE ESTIMATION DEVICE AND SEA STATE ESTIMATION METHOD

Abstract

A sea state estimation device includes a time history memory which stores a time history of ship motion data relating to a rolling angle, a pitching angle and a heaving displacement; a hull condition data calculator which calculates current hull condition data of the hull based on the time history of the ship motion data stored in the time history memory; a cross spectra calculator which calculates cross spectra of respective ship motions; a hull response function calculator which calculates a hull response function based on the hull condition data calculated by the hull condition data calculator; and a sea state estimator which estimates a sea state in sea on which the hull is running based on the cross spectra of the respective ship motions calculated by the cross spectra calculator and the hull response function calculated by the hull response function calculator.

Description

INCORPORATION BY REFERENCE

The present application claims priority from Japanese Patent Application No. 2015-228228 filed on Nov. 20, 2015, the contents of which are hereby incorporated by reference.

TECHNICAL FIELD

The present invention relates to a sea state estimation device and a sea state estimation method.

BACKGROUND ART

Conventionally, with respect to a ship running in irregular waves, it is important to appropriately grasp ship motion data, hull condition data, and a sea state from a safety point of view. The ship motion data are the data related to the motion of the hull, such as a displacement and an acceleration of the hull. Further, the hull condition data are the data related to the condition of the hull, such as a draft, a displacement, and a transverse metacenter height (hereinafter, also referred to as “GM”) of the hull. Further, the sea state is the information related to the sea, such as a wave height, a wave period, and a wave direction in the area where the ship is running.

Conventionally, it is common that the ship motion data of these pieces of information are kinetically obtained based on the previously set hull condition data and the sea state provided from an information providing institution such as the Meteorological Agency.

However, in this method, it is impossible to grasp the ship motion data and the hull condition data appropriately. This is because the sea state provided from an information providing institution such as the Meteorological Agency has a large amount of information on a wide area of ocean instead of a local area of sea in which the ship is currently running, and has a low accuracy.

To address this issue, in recent years, a method is reported in which ship motion data are measured with various devices mounted on a ship, the measured ship motion data, which are unsteady time-series data, are subjected to statistic processing in real time, so that the hull condition data and the sea state are statistically estimated (for example, see Non-Patent Literature 1).

Non-Patent Literature 1 discloses a technique in which the sea state is estimated by analyzing the ship motion data, which are unsteady time-series data, by using a Time Varying Coefficient Vector AR (TVVAR) model.

CITATION LIST

Non Patent Literature

• Non-Patent Literature 1: Tsugikiyo Hirayama, Toshio Iseki, Shigesuke Ishida, “Real-Time Detection Method of Encountering Ocean Waves and Recent Results”, The Japan Society of Naval Architects and Ocean Engineers, December 2003, pp. 74-96.

SUMMARY OF INVENTION

Technical Problems

In the technique disclosed in the above-mentioned Non-Patent Literature 1, a hull response function is used for estimating a sea state. “hull response function” is a function indicating how the hull responds (moves) when the hull receives waves with a regular wavelength from any given direction, and the parameters of the function are a wave direction, a wave length and so on. In the technique disclosed in Non-Patent Literature 1, this hull response function is set based on a ship speed which is a fixed value.

However, a ship speed changes from moment to moment on the actual sea and hence, accuracy of the hull response function decided based on a ship speed which is a fixed value is low. As a result, eventually, accuracy in the estimation of a sea state is also low.

The present invention has been made in view of the above-mentioned circumstances, and it is an object of the present invention to provide a sea state estimation device and a sea state estimation method which can estimate a sea state with high accuracy in response to a change in a ship speed.

Solution to Problems

To achieve the above-mentioned object of the present invention, a sea state estimation device according to the present invention includes: a time history memory unit which stores a time history of ship motion data relating to a rolling angle, a pitching angle and a heaving displacement of a ship; a hull condition data calculation unit which calculates current hull condition data of the ship based on the time history of the ship motion data stored in the time history memory unit; a cross spectra calculation unit which calculates cross spectra of respective ship motions consisting of rolling, pitching and heaving of the ship based on the time history of the ship motion data stored in the time history memory unit; a hull response function calculation unit which calculates a hull response function based on the hull condition data calculated by the hull condition data calculation unit; and a sea state estimation unit which estimates a sea state in sea on which the ship is running based on the cross spectra of the respective ship motions calculated by the cross spectra calculation unit and the hull response function calculated by the hull response function calculation unit.

To achieve the above-mentioned object of the present invention, a sea state estimation method according to the present invention is a sea state estimation method in a sea state estimation device which includes a time history memory unit which stores a time history of ship motion data relating to a rolling angle, a pitching angle and a heaving displacement of a ship, the method including: a hull condition data calculation step where a hull condition data calculation unit calculates current hull condition data of the ship based on the time history of the ship motion data stored in the time history memory unit; a cross spectra calculation step where a cross spectra calculation unit calculates cross spectra of respective ship motions consisting of rolling, pitching and heaving of the ship based on the time history of the ship motion data stored in the time history memory unit; a hull response function calculation step where a hull response function calculation unit calculates a hull response function based on the hull condition data calculated in the hull condition data calculation step; and a sea state estimation step where a sea state estimation unit estimates a sea state in sea on which the ship is running based on the cross spectra of the respective ship motions calculated in the cross spectra calculation step and the hull response function calculated in the hull response function calculation step.

Advantageous Effects of Invention

According to the present invention, a sea state can be estimated with high accuracy corresponding to a change in a ship speed.

BRIEF DESCRIPTION OF DRAWINGS

FIG. 1 is a block diagram showing an example of a hardware configuration of a sea state estimation system including a sea state estimation device according to an embodiment.

FIG. 2 is a block diagram showing a functional configuration example of the sea state estimation system including the sea state estimation device according to the embodiment.

FIG. 3 is a flowchart showing one example of a control logic of the sea state estimation system including the sea state estimation device according to the embodiment

FIG. 4 is a table showing one example of data structure of a data base which a hull response function calculation section according to the embodiment uses.

FIG. 5A is a graph for explaining a hull response function according to the embodiment.

FIG. 5B is a graph for explaining a hull response function according to the embodiment.

FIG. 5C is a graph for explaining a hull response function according to the embodiment.

FIG. 6A is a graph showing one example of ship motion data according to the embodiment.

FIG. 6B is a graph showing one example of ship motion data according to the embodiment.

FIG. 6C is a graph showing one example of ship motion data according to the embodiment.

FIG. 7A is a graph showing one example of cross spectra according to the embodiment.

FIG. 7B is a graph showing one example of cross spectra according to the embodiment.

FIG. 7C is a graph showing one example of cross spectra according to the embodiment.

FIG. 7D is a graph showing one example of cross spectra according to the embodiment.

FIG. 7E is a graph showing one example of cross spectra according to the embodiment.

FIG. 7F is a graph showing one example of cross spectra according to the embodiment.

FIG. 8A is a graph showing one example of a directional wave spectrum according to the embodiment.

FIG. 8B is a graph showing one example of a directional wave spectrum according to the embodiment.

DESCRIPTION OF EMBODIMENT

Hereinafter, an embodiment according to the present invention will be described.

[Hardware Configuration of Sea State Estimation Device]

FIG. 1 is a block diagram showing an example of a hardware configuration of a sea state estimation system including a sea state estimation device according to the embodiment.

A sea state estimation system 1 shown in FIG. 1 is equipped with a satellite compass 2, an information processing device (sea state estimation device) 3, and a display 4. The sea state estimation system 1 is mounted in a hull of a ship (vessel).

The satellite compass (GPS compass) 2 is a device having a function as a direction sensor which calculates the direction of the ship based on the relationship between relative positions of two GPS antennas attached in a bow direction of the ship. The satellite compass 2 has also a function as a motion-in-wave sensor which can measure a transverse motion in wave (rolling), a longitudinal motion in wave (pitching), and a vertical motion in wave (heaving) of the ship. A gyro sensor may be used instead of the satellite compass 2.

The information processing device 3 is a computer device which includes a memory device 31, a processing unit 32, an interface device 33, an input device 34, an auxiliary storage device 35, and a drive device 36 which are each connected to one another through a bus 38. The information processing device 3 estimates the sea state based on information measured by the satellite compass 2. The information processing device 3 corresponds to a “sea state estimation device” of the claims. The information processing device 3 and the display 4 to be described later may be integrated with the satellite compass 2.

The memory device 31 is a storage device such as a random access memory (RAM) which reads out and stores, at a time of start-up of the information processing device 3, a program (a program which realizes functions of the sections ranging from the cross spectra calculation section 23 to the sea state estimation section 26 shown in FIG. 2) and the like stored in the auxiliary storage device 35. The memory device 31 stores also a file, data and the like necessary to execute the program.

The processing unit 32 is a processing unit such as a central processing unit (CPU) which executes a program stored in the memory device 31. The interface device 33 is an interface device to connect the processing unit 32 to an external device such as the satellite compass 2 and the display 4. The input device 34 is an input device (for example, a keyboard or a mouse) to provide a user interface.

The auxiliary storage device 35 is a storage device such as a hard disk drive (HDD) which stores a program, a file, data, and the like. The auxiliary storage device 35 stores a program and the like which realize functions of the sections ranging from the cross spectra calculation section 23 to the sea state estimation section 26 shown in FIG. 2.

The drive device 36 is a device which reads out a program (for example, a program which realizes functions of the sections ranging from the cross spectra calculation section 23 to the sea state estimation section 26 shown in FIG. 2) stored in a storage medium 37. A program read out by the drive device 36 is installed in the auxiliary storage device 35. The storage medium 37 is a storage medium such as a universal serial bus (USE) memory and an SD memory card, and the storage medium 37 stores the above program or the like.

The display 4 is an output device which outputs, on a screen, output data generated by the information processing device 3, for example, the sea state.

[Functional Configuration of System]

FIG. 2 is a block diagram showing a functional configuration example of the sea state estimation system including the sea state estimation device according to the embodiment. In the following description, components similar to the components shown in FIG. 1 are assigned the same reference numerals, and the repeated description is appropriately omitted.

The sea state estimation system 1 shown in FIG. 2 has a measurement section 21, a time history memory 22, a cross spectra calculation section 23, a hull condition data calculation section 24, a hull response function calculation section 25, the sea state calculation section 26, a response function memory 27, and an output section 28. The sea state estimation system 1 is mounted on a ship. The information processing device 3 realizes respective functions of the time history memory 22, the cross spectra calculation section 23, the hull condition data calculation section 24, the hull response function calculation section 25, the sea state calculation section 26, and the response function memory 27.

The measurement section 21 is a measurement unit which measures ship motion data of a ship on which the sea state estimation system 1 is mounted. Here, the ship motion data means the data related to the motion of the ship such as a ship speed, a rolling angle, a pitching angle, a displacement of heaving, and the like of the hull. It is also possible to use angular velocities of rolling and pitching and an acceleration of heaving of the ship. The measurement section 21 is realized by the satellite compass 2 shown in FIG. 1 or a gyro sensor.

The time history memory 22 is a time history memory unit which stores a time history of the ship motion data measured by the measurement section 21. The time history memory 22 stores time-series data of the ship motion data in a predetermined period from past to now. The time history memory 22 is realized by the memory device 31 and the like shown in FIG. 1. An input source of the time history of the ship motion data is not limited to the measurement section 21. For example, it is also possible that another information processing device storing a time history of the ship motion data is used as the input source.

The cross spectra calculation section 23 calculates cross spectra of the respective ship motions (the rolling, the pitching, and the heaving) based on the time-series data of the ship motion data in a predetermined period stored in the time history memory 22. The cross spectra calculation section 23 is realized by the processing unit 32 and the like shown in FIG. 1.

The hull condition data calculation section 24 is a hull condition data calculation unit which calculates the hull condition data based on the time-series data which is the time history of the ship motion data stored in the time history memory 22 in a predetermined period. Here, the hull condition data means the data related to the condition of the hull such as a draft, a displacement, and a GM of the hull. The hull condition data calculation section 24 is realized by the processing unit 32 or the like shown in FIG. 1.

The hull response function calculation section 25 is a hull response function calculation unit which calculates a hull response function of the hull based on the hull condition data calculated by the hull condition data calculation section 24 and the data relating to the hull response function stored in the response function memory 27. The hull response function calculation section 25 is realized by the processing unit 32 and the like shown in FIG. 1.

The sea state estimation section 26 is a sea state estimation unit which estimates a local sea state in the sea in which the ship equipped with the sea state estimation system 1 is running, based on the cross spectra calculated by the cross spectra calculation section 23 and the hull response function calculated by the hull response function calculation section 25. Here, the sea state means the information related to the sea state such as a wave height, a wave period, a wave direction, and the like of the wave in the area in which the ship is running.

The sea state estimation section 26, firstly, probability statistically calculates a directional wave spectrum based on the cross spectra and the hull response function. Then, the sea state estimation section 26 estimates the sea state based on the calculated directional wave spectrum. The detailed process is described later. The sea state estimation section 26 is realized by the processing unit 32 and the like shown in FIG. 1.

The response function memory 27 is a response function storage unit which stores data related to the hull response function which is used when the hull response function calculation section 25 calculates the hull response function. “hull response function” is a function indicating how the ship responds (moves) when the ship receives waves with a regular wavelength from any given direction, and the parameters of the function are a wave direction, a wave length, and the like. “data related to the hull response function” is a data base in which the hull response function of the ship relating to regular waves is made to correspond to the hull condition data (the draft and the GM), a ship speed, and a damping coefficient. The response function memory 27 is realized by the memory device 31 and the like shown in FIG. 1.

The output section 28 is an output unit which outputs the hull condition data calculated by the hull condition data calculation section 24 and the sea state estimated by the sea state estimation section 26. The output section 28 is realized by the display 4 and the like shown in FIG. 1.

With the configuration described above, in the sea state estimation system 1 according to the embodiment, the sea state estimation section 26 estimates the sea state based on the ship motion data measured by the measurement section 21. An output section 28 outputs the estimated sea state.

[Control Logic of System]

FIG. 3 is a flowchart showing a control logic of the sea state estimation system including the sea state estimation device according to the embodiment.

The sea state estimation system 1 repeatedly estimates the sea state by repeatedly performing a control logic of a series of steps S1 to S8 shown in FIG. 3. A description will be given below, appropriately with reference to FIG. 2.

First, in step S1, the measurement section 21 measures the ship motion data (step S1). Specifically, the data of the rolling angle, the pitching angle, and the displacement of heaving of the hull are measured. Step S1 may be successively performed in the process of repeating the series of steps S1 to S8, or may be repeatedly performed, by a batch process or the like, as a process independent from steps S2 to S8.

By repeatedly performing the process of step S1 shown in FIG. 3, the time-series data of the rolling angle, the pitching angle, and the displacement of heaving in a predetermined period from past to now, in other words, the data of the angles and the like at respective times are stored in the time history memory 22. For example, ship motion data is measured and stored at an interval of 1 second to 10 seconds, and respective data are expressed by a value which includes the second decimal place or more decimal place. “ship motion data” includes unsteady time-series data which includes an effect exerted by a change in course or the acceleration or deceleration of a speed besides motions due to a sea condition state. In step S1, such data are separated and are not used in succeeding steps. The Akaike's Information Criterion (AIC) method which is an existing technique is used in this process. Then, the process goes to respective processes of steps S2, S4, and S6.

Next, in step S2, the hull condition data calculation section 24 calculates (estimates) a rolling natural frequency based on the time-series data of the rolling angle (step S2). Although the process of step S2 is a known technique, an example of the process will be described below.

Specifically, a second order linear probability dynamic model (see Equation (1) below) about time-series data x(t) of rolling is considered. In Equation (1), a₁ (=2α) is a damping coefficient and a₂ (=ω²) is a square of a natural angular frequency ω. The term u(t) represents an external force term dealt as a stochastic process and has a finite dispersion. However, the term u(t) does not need whiteness.

[Mathematical Expression 1]

x″(t)+a₁x′(t)+a₂x(t)=u(t)  (1)

Further, the external force term u(t) in Equation (1) is expressed by an m order continuous auto regression model shown in Equation (2) below. In Equation (2), bi (i=1, . . . , m) is a coefficient of the model, and v(t) is a normal white noise, where the average is 0, and the dispersion is σ².

$\begin{array}{cc}\left[\mathrm{Mathematical}\mathrm{Expression}2\right]& \\ u^{\left(m\right)}\left(t\right)+\sum _{i=1}^mb_iu^{\left(m-i\right)}\left(t\right)=v\left(t\right)& \left(2\right)\end{array}$

By substituting Equation (1) into Equation (2), a whitened (m+2) order continuous auto regression model shown in Equation (3) is obtained. c_i (i=1, . . . , m+2) in Equation (3) is a coefficient of the model.

$\begin{array}{cc}\left[\mathrm{Mathematical}\mathrm{Expression}3\right]& \\ u^{\left(m+2\right)}\left(t\right)+\sum _{i=1}^{m+2}c_iu^{\left(m+2-i\right)}\left(t\right)=v\left(t\right)& \left(3\right)\end{array}$

Equation (3) is expressed as Equation (4) below in a vector form.

$\begin{array}{cc}\left[\mathrm{Mathematical}\mathrm{Expression}4\right]& \\ x^{\prime }\left(t\right)=\mathrm{Ax}\left(t\right)+\mathrm{Bv}\left(t\right)\mathrm{where}\{\begin{array}{c}x\left(t\right)={\left(x\left(t\right)x^{\prime }\left(t\right)x^{″}\left(t\right)\dots x^{\left(m+1\right)}\right)}^T\\ A=\left(\begin{array}{ccccc}0& 1& 0& \dots & 0\\ 0& 0& 1& 0& ⋮\\ ⋮& ⋮& \ddots & \ddots & 0\\ 0& 0& \dots & 0& 1\\ -C_{m+2}& -C_{m+1}& \dots & -C_2& -C_1\end{array}\right)\\ B=I_{\left(m+2\right)\times \left(m+2\right)}\\ v\left(t\right)={\left(00\dots 0v\left(t\right)\right)}^T\end{array}& \left(4\right)\end{array}$

Equation (4) is dealt as a system model of a state space model after being discretized. Further, in order to simultaneously estimate the unknown coefficients c_i (where i=1, . . . , m+2), a state vector of the state space model is considered with the above-described unknown coefficients a_i, b_i included therein, and the model is then extended into an auto-organizing state space model; then, the state estimation and the unknown coefficients are simultaneously estimated by using an Ensemble Kalman Filter. The state estimation by using the Ensemble Kalman Filter is a known technique and is thus not described here.

By the procedure described above, in step S2, the current rolling natural frequency is calculated (estimated) based on the time-series data of the rolling angle stored in the time history memory 22 in a predetermined period. The rolling natural frequency may be calculated by a method different from the above method. For example, the rolling natural frequency may be calculated by using a discretized auto regression model.

After that, in step S3, the hull condition data calculation section 24 calculates (estimates) the GM based on the rolling natural frequency calculated in step S2 (step S3).

In step S3, in a non-linear observation model in which the rolling natural frequency (or a rolling natural period, which is the inverse) calculated in step S2 is used as observed data and in which the GM and the radius of gyration are the state variables, it is assumed that the state variables fluctuate slightly with time, and this is considered as the system model; thus, a general state-space model analysis is performed to simultaneously estimate the GM and the radius of gyration.

That is, the relationship represented by Equation (5) below holds between the rolling natural frequency f (=ω/2π) calculated in step S2 and the GM. In Equation (5), T is the rolling natural period, f is the rolling natural frequency, k is the radius of gyration, and g is the gravitational acceleration. In Equation (5), GM and k are unknown.

$\begin{array}{cc}\left[\mathrm{Mathematical}\mathrm{Expression}5\right]& \\ T=\frac{1}{f}=\frac{2\pi k}{\sqrt{\mathrm{gGM}}}& \left(5\right)\end{array}$

Then, in step S3, the state space model is considered in which two unknowns of the GM and the radius of gyration k are the state variables. That is, a non-linear observation model is considered in which the rolling natural period T_n is observed from the estimated amount of the states of GM_n and k_n at time n. It is assumed that the state variables fluctuate slightly with time, and this is considered as the system model. Thus, the general state-space model represented by Equation (6) below is configured. In Equation (6), v_n is a system noise, and w_n is an observation noise. For the sake of simplicity, the noises are considered as normal white noises.

$\begin{array}{cc}\left[\mathrm{Mathematical}\mathrm{Expression}6\right]& \\ \{\begin{array}{c}{x}_n=x_{n-1}+v_n\\ T_n=h\left(x_n,w_n\right)\end{array}\mathrm{where}\{\begin{array}{c}{x}_n={\left({\mathrm{GM}}_nk_n\right)}^T\\ v_n={\left(v_{1,n}v_{2,n}\right)}^T\end{array}& \left(6\right)\end{array}$

In step S3, state estimation is performed by using a Monte Carlo filter, which is a type of particle filters, based on the general state-space model represented by Equation (6). The state estimation by using a Monte Carlo filter is a known technique and is not described here.

By the process of steps S2 and S3 described above, the hull condition data calculation section 24 calculates (estimates) the current GM of the ship in real time by analyzing the time-series data of the rolling angle in a predetermined period stored in the time history memory 22. By this method, it is possible to grasp the change in the position of the center of gravity in a loading state and grasp a degree of motion of the ship in real time.

Further, because no approximation is used to estimate GM, it is possible to obtain a generally stable estimation result of GM. For example, in the case that GM of a ship at the time of design was 0.52, GM was estimated to be in the range of 0.48 to 0.54 by the method according to the embodiment. Therefore, the method according to the embodiment enables GM to be estimated with high accuracy. In the above Equation (6), the equation may be made to include a variable older than the previous timing (for example X_n-2).

Returning to FIG. 3, if the process goes from step S1 to step S4, the hull condition data calculation section 24 calculates the current draft and the current displacement of the ship by analyzing the time-series data of the displacement of heaving stored in the time history memory 22 in step S1 (step S4).

In step S4, the hull condition data calculation section 24 calculates (estimates) a draft based on tidal data which is calculated based on, for example, an installation height of a GPS antenna (corresponding to the measurement section 21 shown in FIG. 1), an inclination angle of the ship in the longitudinal direction, and Table on tidal harmonic constant along coasts of Japan published by Japan Coast Guard. Further, a natural frequency of heaving is calculated (estimated) by discretizing a motion equation of heaving prepared based on time-series data of displacement of heaving in the same manner as rolling in step S1. Then, a displacement is calculated (estimated) by making use of a relationship between the estimated natural frequency, a TPC (Tons Per Centimeter) which corresponds to the combination of a draft preliminarily calculated based on a hull diagram and an inclination angle of the ship in the longitudinal direction and a displacement.

Further, in step S4, for example, a displacement and a trim are calculated based on time-series data on a pitching angle and displacement of heaving stored in the time history memory 22 using a statistical analysis method. With such a calculation, drafts at a bow and a stern can be decided using a hydro calculation table of the hull. When the ship in which the sea state estimation system 1 according to the embodiment is mounted includes a draft meter, a displacement is calculated based on data measured using the draft meter. When the ship does not include a draft meter, in the same manner as steps S1 and S2, a displacement is calculated using an auto-organizing state space model analyzing method.

If the process goes from steps S3 and S4 to step S5, the hull response function calculation section 25 calculates a current hull response function of the ship based on the hull condition data calculated in step S3 and S4 (step S5).

As the process of step S5, there are two possible methods described below, and any one of them is used. In the first method, a hull response function is previously calculated to make a data base (for example, FIG. 4) by using as parameters the hull condition data (the draft and GM), a ship speed, and a sea state (a wave height, a wave period, and a wave direction) which is to be an input; and the hull response function corresponding to the current states is then obtained by an interpolation calculation such as a linear interpolation, a quadratic interpolation or the like. In the second method, the hull response function is directly calculated based on a calculation formula, with the current state of the hull and the sea state being used as inputs. In any of the methods, the sea state is an unknown to be obtained below and is also a term necessary to calculate the hull response function. The most appropriate response function is selected in real time as a non-linear problem by using an iteration method. In this specification, as shown in FIG. 5, “hull response function” is a result of calculation of waves which come from respective directions (from 0 degree (following) to 180 degrees (ahead) at an arbitrary wave height, that is, the expression of a characteristic of a hull response in the form of function, and differs from a hull condition. “hull response” expresses an effect which the ship receives from an external force such as a wave, and is expressed as a ship motion ora change in structure (structural response).

FIG. 4 is a Table showing one example of data structure of a data base which a hull response function calculation section according to the embodiment uses.

In the data base 41 shown in FIG. 4, GM (1≤k≤o), a displacement (1≤j≤n) and hull response functions f_xyz (1≤x≤o, 1≤y≤n, 1≤z≤m) for respective ship speeds (1≤g≤m) are stored. Index h, p, r in the hull response function f_xyz indicate heaving, pitching and rolling respectively. Such data base is stored in a response function memory 27 (see FIG. 2). In step S5, the hull response function calculation section 25 obtains a hull response function which corresponds to hull condition data calculated in steps S3 and S4 by selecting from the data base or by selecting from the data base and by performing an interpolation such as a linear interpolation or a quadratic interpolation.

FIG. 5A to FIG. 5C are figures for explaining a hull response function according to the embodiment.

In FIG. 5A, a motion amplitude of heaving (Amp.)/(wave height/2)(Aw) is taken on an axis of ordinates and a wavelength/length between perpendiculars (Lpp) of heaving is taken on an axis of abscissas. In FIG. 5B, a motion amplitude of pitching (Amp.)/(wave height/2)(Aw) is taken on an axis of ordinates and a wavelength/length between perpendiculars (Lpp) of pitching is taken on an axis of abscissas. In FIG. 5C, a motion amplitude of rolling (Amp.)/(wave height/2)(Aw) is taken on an axis of ordinates and a wavelength/length between perpendiculars (Lpp) of rolling is taken on an axis of abscissas. In FIG. 5A to FIG. 5C, a calculation example of a motion amplitude is shown with respect to the response functions f_xyz of heaving, pitching and rolling of the hull respectively while using a sea state, a wave height, a wavelength (λ), and a wave direction (ω) as parameters.

In the process of step S5, the most appropriate hull response function is selected in real time, depending on the current ship motions (the pitching, the rolling, and the heaving) and the ship speed measured by the measurement section 21 (the satellite compass 2). By this process, even on the actual sea in which the ship speed changes from moment to moment, the most appropriate response function of the hull can be obtained. As a result, it is also possible to improve the accuracy of the estimation of the sea state to be described later.

Further, if the process goes from step S1 to step S6, the cross spectra calculation section 23 calculates cross spectra of the respective ship motions (the rolling, the pitching, and the heaving) based on the time-series data of the displacement of heaving, the pitching angle, and the rolling angle stored in the time history memory 22 in step S1 (step S6).

In step S6, based on the time-series data of the displacement of heaving (unit: m), the pitching angle (unit: rad), the rolling angle (unit: rad), the cross spectra calculation section 23 calculates for each frequency, cross spectra of the respective ship motions (the rolling, the pitching, and the heaving) including a rolling auto spectrum (unit: rad²/s), a pitching auto spectrum (unit: rad²/s), a heaving auto spectrum (unit: m²/s), a pitching-heaving cross spectrum (unit: rad·m/s), a rolling-heaving cross spectrum (unit: rad·m/s), and a pitching-rolling cross spectrum (unit: rad²/s). The cross spectra obtained for respective frequencies are stored as the time-series data in the time history memory 22.

In step S6, the cross spectra calculation section 23 calculates cross spectra using an auto regression model analyzing method. That is, an auto regression model is applied to time-series data of displacement of heaving, a pitching angle and rolling angle, and cross spectra are calculated using an estimated auto regression coefficient and a variance-covariance matrix. The determination of model order is important in such a calculation, and the model order is determined automatically using an AIC method. These series of procedures are referred to as a MAICE (Minimum AIC Estimation) method which was established in 1980s.

If the process goes from steps S5 and S6 to step S7, the sea state estimation section 26 probability statistically calculates a directional wave spectrum based on the current hull response functions calculated in step S5 and the cross spectra of the respective ship motions (the rolling, the pitching, and the heaving) calculated in step S6 (step S7).

FIG. 6A to FIG. 6C are figures showing one example of ship motion data according to the embodiment.

FIG. 6A to FIG. 6C respectively show ship motion data in waves measured in step S1 with respect to a ship having a length of approximately 70 m. That is, FIG. 6A to FIG. 6C respectively show time-series data of a pitching angle, a rolling angle and the displacement of heaving.

FIG. 7A to FIG. 7F are graphs showing one example of cross spectra according to the embodiment.

FIG. 7A to FIG. 7F show one example of a result of an analysis of cross spectra calculated in step S6. That is, FIG. 7A to FIG. 7F show a pitching auto spectrum, pitching-rolling cross spectra, pitching-heaving cross spectra, a rolling auto spectrum, rolling-heaving cross spectra, and a heaving auto spectrum.

FIG. 8A and FIG. 8B are figures showing one example of the directional wave spectrum according to the embodiment.

FIG. 8A and FIG. 8B respectively show a one-dimensional spectrum and a direction distribution of a directional wave spectrum obtained by reverse calculation in step S7. A sea state, a wave height, a wave period and a wave direction are obtained based on this one-dimensional spectrum and the direction distribution. Hereinafter, such process will be described together with a theory according to the embodiment used in step S7.

If it is assumed that an ocean wave is represented by superposing component waves coming from every direction and having all frequencies, a sea surface variation amount η(t) at time t at an fixed point (ship position) is expressed by Equation (7) below using a directional wave spectrum E(f, x) (unit: m²/(rad/s)). In Equation (7), the part under the root symbol and ε(f, x) are respectively the amplitude and the phase of a component wave having a frequency f coming from a direction x.

[Mathematical Expression 7]

n(t)=∫_−π^π∫₀^∞ cos {2πft+ε(f,x)}√{square root over (2E(f,x)dfdx)}  (7)

On the other hand, if it is assumed that the ship motion in wave responds linearly to an input wave, a relationship between the directional wave spectrum E(f_e, x) at an encounter frequency f_e of one wave, and a cross spectrum ϕ_ln(f_e) of the ship motion in wave is expressed generally by Equation (8) below. In Equation (8), l and n are the modes of the ship motion in wave, and H₁(f_e, x) and H_n*(f_e, x) are respectively the response functions in the modes l and n. Further, x is an encounter angle with respect to a wave, and the symbol “*” represents a complex conjugate.

[Mathematical Expression 8]

ϕ_ln(f_e)=∫_−π^πH₁(f_e,x)H*_n(f_e,x)E(f_e,x)dx  (8)

Because Equation (8) is expressed based on the encounter frequency, this equation is converted into an equation based on an absolute frequency (see Equation (9) below).

$\begin{array}{cc}\left[\mathrm{Mathematical}\mathrm{Expression}9\right]& \\ {\phi }_{\ln }\left(f_e\right)={\int }_{\frac{n}{2}}^{\pi }H_l\left(f_{01},x\right)H_n^{*}\left(f_{01},x\right)E\left(f_{01},x\right)\frac{{\mathrm{df}}_{01}}{{\mathrm{df}}_e}\mathrm{dx}+{\int }_{-\frac{\pi }{2}}^{\frac{\pi }{2}}H_l\left(f_{01},x\right)H_n^{*}\left(f_{01},x\right)E\left(f_{01},x\right)\frac{{\mathrm{df}}_{01}}{{\mathrm{df}}_e}\mathrm{dx}\left(f_e<\frac{1}{4A}\right)+{\int }_{-\frac{\pi }{2}}^{\frac{\pi }{2}}H_l\left(f_{02},x\right)H_n^{*}\left(f_{02},x\right)E\left(f_{02},x\right)\frac{{\mathrm{df}}_{02}}{{\mathrm{df}}_e}\mathrm{dx}\left(f_e<\frac{1}{4A}\right)+{\int }_{-\frac{\pi }{2}}^{\frac{\pi }{2}}H_l\left(f_{03},x\right)H_n^{*}\left(f_{03},x\right)E\left(f_{03},x\right)\frac{{\mathrm{df}}_{03}}{{\mathrm{df}}_e}\mathrm{dx}+{\int }_{-\pi }^{-\frac{\pi }{2}}H_l\left(f_{01},x\right)H_n^{*}\left(f_{01},x\right)E\left(f_{01},x\right)\frac{{\mathrm{df}}_{01}}{{\mathrm{df}}_e}\mathrm{dx}& \left(9\right)\end{array}$

In Equation (9), the second to fourth terms on the right side represent contribution at a time of a following wave, in other words, represent the degree of the frequency component of the wave when running on the following wave included in the cross spectrum. The parameter A, the three encounter frequencies f₀₁, f₀₂, and f₀₃ corresponding to the absolute frequency, and the Yacoubian are each defined as shown in Equation (10) below. In Equation (10), U is the ship speed and g is the gravitational acceleration.

$\begin{array}{cc}\left[\mathrm{Mathematical}\mathrm{Expression}10\right]& \\ \{\begin{array}{c}A=\frac{2\pi }{g}\bigcup \cos x\\ f_{01}=\frac{1-\sqrt{1-4Af_e}}{2A}\\ f_{02}=\frac{1+\sqrt{1-4Af_e}}{2A}\\ f_{03}=\frac{1+\sqrt{1+4Af_e}}{2A}\\ \frac{{\mathrm{df}}_{01}}{{\mathrm{df}}_e}=\frac{1}{\sqrt{1-4Af_e}}\\ \frac{{\mathrm{df}}_{02}}{{\mathrm{df}}_e}=\frac{1}{\sqrt{1-4Af_e}}\\ \frac{{\mathrm{df}}_{03}}{{\mathrm{df}}_e}=\frac{1}{\sqrt{1+4Af_e}}\end{array}& \left(10\right)\end{array}$

Here, in the case that the integration range with respect to the encounter angle x is divided into a sufficiently large number K of fine interval, the response function and the directional wave spectrum of a variation amount can be constant in each of the fine integral interval. Therefore, Equation (9) can be discretized into Equation (11) below. In Equation (11), K1 (where, 0≤K1≤K/2) is the number of the fine interval which are in a following wave state in the discrete integral range.

$\begin{array}{cc}\left[\mathrm{Mathematical}\mathrm{Expression}11\right]& \\ {\phi }_{ln}\left(f_e\right)=\Delta x\sum _{k=1}^KH_{l,k}\left(f_{01}\right)H_{n,k}^{*}\left(f_{01}\right)E_k\left(f_{01}\right)+\Delta x\sum _{k=1}^{K_1}H_{l,k}\left(f_{02}\right)H_{n,k}^{*}\left(f_{02}\right)E_k\left(f_{02}\right)+\Delta x\sum _{k=1}^{K_1}H_{l,k}\left(f_{03}\right)H_{n,k}^{*}\left(f_{03}\right)E_k\left(f_{03}\right)\mathrm{where}\{\begin{array}{c}\Delta x=\frac{2\pi }{K}\\ E_k\left(f_{0*}\right)=E_k\left(f_{0*},x_k\right)\\ x_k=-\pi +\left(k-1\right)\Delta x\\ H_{l,k}\left(f_{0*}\right)=H_l\left(f_{0*},x_k\right)\\ H_{n,k}^{*}\left(f_{0*}\right)=H_n^{*}\left(f_{0*},x_k\right)\end{array}& \left(11\right)\end{array}$

Here, in the case that the pitching angle, the rolling angle, and the heaving displacement are respectively any given variation amounts θ, ϕ, and η, the cross spectrum Φ(f_e) is a 3×3 matrix, and Equation (11) can be expressed in a matrix as Equation (12) below. In Equation (12), H(f₀₁) is a 3×K matrix, H(f₀₂) and H(f₀₃) are 3×K1 matrices, E(f₀₁) is K×K diagonal matrix, E(f₀₂) and E(f₀₃) are K1×K1 diagonal matrices. Further, the symbol T represents a transposed matrix.

$\begin{array}{cc}\left[\mathrm{Mathematical}\mathrm{Expression}12\right]& \\ \Phi \left(f_e\right)=H\left(f_{01}\right)E\left(f_{01}\right){H\left(f_{01}\right)}^{*T}+H\left(f_{02}\right)E\left(f_{02}\right){H\left(f_{02}\right)}^{*T}+H\left(f_{03}\right)E\left(f_{03}\right){H\left(f_{03}\right)}^{*T}\mathrm{where}\Phi \left(f_e\right)=\left(\begin{array}{ccc}{\Phi }_{\mathrm{\theta \theta }}\left(f_e\right)& {\Phi }_{\mathrm{\theta \phi }}\left(f_e\right)& {\Phi }_{\mathrm{\theta \eta }}\left(f_e\right)\\ {\Phi }_{\mathrm{\phi \theta }}\left(f_e\right)& {\Phi }_{\mathrm{\phi \phi }}\left(f_e\right)& {\Phi }_{\mathrm{\phi \eta }}\left(f_e\right)\\ {\Phi }_{\mathrm{\eta \theta }}\left(f_e\right)& {\Phi }_{\mathrm{\eta \phi }}\left(f_e\right)& {\Phi }_{\mathrm{\eta \eta }}\left(f_e\right)\end{array}\right)H\left(f_{01}\right)=\left(\begin{array}{ccc}{H}_{\theta 1}\left(f_{01}\right)& \dots & H_{\theta K}\left(f_{01}\right)\\ H_{\phi 1}\left(f_{01}\right)& \dots & H_{\phi K}\left(f_{01}\right)\\ H_{\eta 1}\left(f_{01}\right)& \dots & H_{\eta K}\left(f_{01}\right)\end{array}\right),E\left(f_{01}\right)=\left(\begin{array}{ccc}{E}_1\left(f_{01}\right)& \dots & 0\\ ⋮& \ddots & ⋮\\ 0& \dots & E_K\left(f_{01}\right)\end{array}\right)H\left(f_{0i}\right)=\left(\begin{array}{ccc}{H}_{\theta 1}\left(f_{0i}\right)& \dots & H_{\theta K1}\left(f_{0i}\right)\\ H_{\phi 1}\left(f_{0i}\right)& \dots & H_{\phi K1}\left(f_{0i}\right)\\ H_{\eta 1}\left(f_{0i}\right)& \dots & H_{\eta K1}\left(f_{0i}\right)\end{array}\right),E\left(f_{0i}\right)=\left(\begin{array}{ccc}{E}_1\left(f_{0i}\right)& \dots & 0\\ ⋮& \ddots & ⋮\\ 0& \dots & E_{K1}\left(f_{0i}\right)\end{array}\right)i=2,3& \left(12\right)\end{array}$

Because the cross spectrum matrix Φ(f_e) is an Hermitian matrix, it is enough to deal with the upper triangular matrix. Further, in the case that Equation (12) is expressed by separating the real part and the imaginary part and by introducing an error term W associated with observation, Equation (12) can be expressed by a linear regression model represented by Equation (13) below.

[Mathematical Expression 13]

y=Ax+W  (23)

In Equation (13), y is a vector constituted by the real part and the imaginary part of the cross spectrum matrix Φ(f_e). The parameter A is a coefficient matrix constituted by a logical value of the response function of the ship motion in wave. The error term W is a white noise having statistical characteristics of the average 0 and following a variance-covariance matrix Σ. The encounter angle x is an unknown vector constituted by the discretized directional wave spectrum.

In Equation (13), it is assumed that the cross spectrum is obtained in time series, the directional wave spectrum can be estimated in time series based on such an assumption. This corresponds to considering Equation (13) as a time varying system, and Equation (13) can thus be extended into Equation (14) below, where time is represented by suffix t.

[Mathematical Expression 14]

y _t =A _t x _t +W _t  (14)

Equation (14) is formally equivalent to an observation model in a general state-space model. Therefore, by introducing as a system model a smoothing prior distribution, in which the directional wave spectrum changes smoothly with time, (see Equation (15) below), the problem of estimating the directional wave spectrum can be treated as a problem of the state estimation of the general state-space model shown by Equation (15) below.

$\begin{array}{cc}\left[\mathrm{Mathematical}\mathrm{Expression}15\right]& \\ \{\begin{array}{c}{x}_t=x_{t-1}+v_t\\ y_t=A_tx_t+W_t\end{array}& \left(15\right)\end{array}$

In Equation (15), x_t is a state vector, v_t is a system noise vector, y_t is an observation vector, A_t is a state transition matrix, and W_t is an observation noise vector. Here, considering that the directional wave spectrum is not negative, the logarithm of the state vector x_t is replaced anew by x_t to deform Equation (15) into a general state-space model represented by Equation (16) below.

$\begin{array}{cc}\left[\mathrm{Mathematical}\mathrm{Expression}16\right]& \\ \{\begin{array}{c}{x}_t=x_{t-1}+v_t\\ y_t=A_tF\left(x_t\right)+W_t\end{array}& \left(16\right)\end{array}$

Here, F(x_t) means that F(x_t) is exponential to all the elements. Further, the elements of the state vector are configured as Equation (17) below.

$\begin{array}{cc}\left[\mathrm{Mathematical}\mathrm{Expression}17\right]& \\ \{\begin{array}{c}F\left(x_t\right)=\exp \left[x_t\right]\\ x_t^T=\left[\ln \left(x_1,t\right),\dots ,\ln \left(x_j,t\right)\right]\\ {F\left(x_t\right)}^T=\exp \left[\ln \left(x_1,t\right),\dots ,\ln \left(x_j,t\right)\right]\\ x_{j,t}=E_{k,t}\left(f_{0i}\right),i=1\sim m,j=1\sim j,j=m\times k,k=1\sim K.\end{array}& \left(17\right)\end{array}$

In Equation (17), m is the number of division of the absolute frequency of a wave. Equation (16) is a non-linear observation model, in other words, a non-linear state space model. Therefore, in order to estimate the state, it is necessary to use a method effective in non-linear filtering. Conventionally, a particle filter is used, but this filter has a very large calculation load. To address this issue, a state estimation method by using an Ensemble Kalman Filter is introduced in the embodiment. However, an Ensemble Kalman Filter cannot be applied to Equation (16) in a non-linear observation model as it is. To solve this problem, the extended state vector represented by Equation (18) below is considered.

$\begin{array}{cc}\left[\mathrm{Mathematical}\mathrm{Expression}18\right]& \\ Z_t=\left(\begin{array}{c}{x}_t\\ A_tF\left(x_t\right)\end{array}\right)& \left(18\right)\end{array}$

Further, the extended observation matrix and the extended state transition vector represented by Equation (19) below are considered.

$\begin{array}{cc}\left[\mathrm{Mathematical}\mathrm{Expression}19\right]& \\ \{\begin{array}{c}{\tilde{A}}_t=\left(O_{l\times k}I_{t\times t}\right),\\ {\stackrel{\_}{f}}_t\left(z_{t-1},v_t\right)=\left(\begin{array}{c}{x}_{t-1}+v_t\\ A_tF\left(x_{t-1}+v_t\right)\end{array}\right)\end{array}& \left(19\right)\end{array}$

As a result, Equation (20) below holds for x_t, and an extended system model is obtained. Further, regarding to y_t, a formally linear extended observation model represented by Equation (21) below can be obtained. Because x_t and y_t are an extended state space model in the linear observation, the state estimation by an Ensemble Kalman Filter can be realized. The application of the Ensemble Kalman Filter is a known technique and is thus not described here.

[Mathematical Expression 20]

z _t ={tilde over (f)} _t(z_t-1,v_t)  (20)

[Mathematical Expression 21]

y _t =Ā _t z _t +W _t  (21)

By the above-described process of step S7, the sea state estimation section 26 probability statistically calculates the directional wave spectrum based on the hull response function calculated in step S5 and the cross spectra of the respective ship motions calculated in step S6 (step S7).

In the process of step S7, by probability statistically processing the hull response function in a predetermined period from past to now and the time-series data of the cross spectra of the respective ship motions, the current directional wave spectrum is calculated in real time. Thus, the directional wave spectrum with high accuracy can be derived.

Further, by this method, the directional wave spectrum is estimated based on the state estimation by using the Ensemble Kalman Filter. Thus, it is possible to estimate the directional wave spectrum with high accuracy in much shorter calculation time than in the method using the conventional Monte Carlo filter.

When the process of step S7 has finished, the process goes to step S8, and the sea state estimation section 26 estimates the sea state based on the directional wave spectrum calculated in step S7 (step S8). In step S8, it is possible to estimate the sea state such as the wave direction, the wave period, the significant wave height, and the like in the local sea in which the ship is running, based on the directional wave spectrum calculated in step S7.

To be more specific, a wave direction is estimated by integrating directional wave spectra with respect to a frequency, and wave periods (an average wave period, a zero-cross wave period, and a wave period between extreme values), and a significant wave height are estimated by integrating directional wave spectra with respect to a directional angle. To acquire a stable and normal result of analysis, a method is used where patterns of spectra of a normal result of analysis and a pattern of spectra of an abnormal result of analysis are recognized and the abnormal result of analysis is screened. As described previously, the calculation of a cross spectrum in step S6 is performed by applying an auto regression model to the calculation. In this case, an AIC is used for determining a model order and hence, stationarity of time sequence can be determined by comparing AIC values (this method being referred to as a local stationary auto regression model analyzing method). Screening of an abnormal result of analysis is performed in accordance with the following steps. That is, in first step, stationarity of time-series data is determined using an AIC. When the time-series data is stationary, the procedure goes to next step (second step). In second step, across spectrum is calculated. Then, in third step, a directional wave spectrum is calculated. Then, in fourth step, a response spectrum is calculated based on the estimated directional wave spectrum and a response function of a motion in wave. Then, in fifth step, the response spectrum calculated in fourth step and an auto component of the cross spectrum calculated in second step are compared to each other. Then, in sixth step, reliability of the estimation result is evaluated in accordance with machine learning based on a differential acquired in fifth step.

An embodiment of the present invention is described above, but the above embodiment describes merely one of the application examples of the present invention. Therefore, it is not intended that the technical scope of the present invention is limited to the specific configuration of the above embodiment.

REFERENCE SIGNS LIST

• • 1: sea state estimation system
  • 2: satellite compass
  • 3: information processing device (sea state estimation device)
  • 4: display
  • 21: measurement section
  • 22: time history memory
  • 23: cross spectra calculation section
  • 24: hull condition data calculation section
  • 25: hull response function calculation section
  • 26: sea state estimation section
  • 27: response function memory
  • 28: output section

Claims

1. A sea state estimation device comprising: a time history memory unit which stores a time history of ship motion data relating to a rolling angle, a pitching angle and a heaving displacement of a hull;a hull condition data calculation unit which calculates current hull condition data of the ship based on the time history of the ship motion data stored in the time history memory unit;a cross spectra calculation unit which calculates cross spectra of respective ship motions consisting of rolling, pitching and heaving of the hull based on the time history of the ship motion data stored in the time history memory unit;a hull response function calculation unit which calculates a hull response function based on the hull condition data calculated by the hull condition data calculation unit; anda sea state estimation unit which estimates a sea state in sea on which the ship is running based on the cross spectra of the respective ship motions calculated by the cross spectra calculation unit and the hull response function calculated by the hull response function calculation unit.

2. The sea state estimation device according to claim 1, wherein the hull condition data includes a draft, a displacement and a transverse metacenter height of the hull.

3. The sea state estimation device according to claim 1 further comprising: a satellite compass; anda response function storage unit which stores a response function of the hull relating to regular waves, and the hull condition data and a ship speed in a sea state where the response function is made to correspond to the hull condition data and the ship speed, andthe hull response function calculation unit selects the response function which is made to correspond to the current hull condition data and the ship speed measured by the satellite compass, and calculates the hull response function of the hull based on the selected response function.

4. The sea state estimation device according to claim 1, wherein the sea state estimation unit calculates probability statistically a directional wave spectrum based on the respective cross spectra of the ship motions and the hull response function, and estimates a sea state based on the calculated directional wave spectrum.

5. A sea state estimation method in a sea state estimation device which includes a time history memory unit which stores a time history of ship motion data relating to a rolling angle, a pitching angle and displacement of heaving of a hull, the method comprising: a hull condition data calculation step where a hull condition data calculation unit calculates current hull condition data of the hull based on the time history of the ship motion data stored in the time history memory unit;a cross spectra calculation step where a cross spectra calculation unit calculates cross spectra of respective ship motions consisting of rolling, pitching and heaving of the hull based on the time history of the ship motion data stored in the time history memory unit;a hull response function calculation step where a hull response function calculation unit calculates a hull response function based on the hull condition data calculated in the hull condition data calculation step; anda sea state estimation step where a sea state estimation unit estimates a sea state in sea on which the hull is running based on the cross spectra of the respective ship motions calculated in the cross spectra calculation step and the hull response function calculated in the hull response function calculation step.