In order to solve the weak nonlinear problem in the simulation of strong nonlinear freak waves, an improved phase modulation method is proposed based on the Longuet-Higgins model and the comparative experiments of wave spectrum in this paper. Experiments show that this method can simulate the freak waves at fixed time and fixed space coordinates. In addition, by comparing the target wave spectrum and the freak wave measured in Tokai of Japan from the perspective of B-F instability and spectral peakedness, it is proved that the waveform of the simulated freak waves can not only maintain the spectral structure of the target ocean wave spectrum, but also accord with the statistical characteristics of the wave sequences. Then, based on the Kirchhoff approximation method and the modified Two-Scale Method, the electromagnetic scattering model of the simulated freak waves is established, and the normalized radar cross section (NRCS) of the freak waves and their background sea surfaces is analyzed. The calculation results show that the NRCS of the freak waves is usually smaller than their large-scale background sea surfaces. It can be concluded that when the neighborhood NRCS difference is less than or equal to −12 dB, we can determine where the freak waves are.

1. Introduction

Freak waves are high and steep waves in the ocean. The duration of freak waves is very short, but the contingency and great destructiveness are extremely threatening to shipping and marine engineering structures. Therefore, the study of freak waves has attracted more and more attention [1]. The occurrence mechanism and engineering prediction of freak waves have become a hot research topic in the field of physical oceanography and ship hydrodynamics [2]. However, the reasons for the occurrence of the freak waves are still unclear, their occurrence has many uncertainties, and they are difficult to observe by the fixed point marine buoys, shore based radars, or optical sensors. It has been confirmed that the formation mechanism of freak waves is divided into linear mechanism and nonlinear mechanism, but the specific reasons are still being explored. Freak waves occur suddenly, have a short duration, and have a wide range of time and space in the global ocean. Therefore, they are very difficult to record. Because most phenomena of the freak waves cannot be observed, which means the data cannot be collected actually, it is very important to use the technology of mobile radar and satellite remote sensing to study and simulate the freak waves. Synthetic Aperture Radar (SAR) is capable of capturing high-resolution microwave images. The microwaves are highly penetrative and can work under any weather conditions. For the safety of maritime navigation and offshore platforms, it is of great significance to explore the physical mechanism of the occurrence, evolution, and extinction of freak waves [36]. Moreover, the use of SAR to monitor and predict freak waves is a disaster reduction technology that should be valued [6].

Given the difficulty of using SAR to observe the freak waves, more and more scientists are beginning to pay attention to numerical simulation methods [7, 8]. Now, numerical simulations of freak waves in deep water are mainly based on linear superposition methods or cubic nonlinear Schrodinger equation. Based on nonlinear wave equation of wave modulation instability, we can study the occurrence mechanism of freak waves [911]. However, due to the large amount of computation, this method is not easy to apply in engineering. At the same time, it is very difficult to control the time and space conditions in the simulation of freak waves [1215]. The method based on the Longuet-Higgins model is effective in simulating the freak waves in the laboratory, which is simple and practical [16]. Lawton simulated the freak waves to form the random initial phase of waves through artificial intervention, which was an inefficient method, and we could not control the generation time and place of freak waves [17]. Pei simulated the recorded freak waves using the three wave trains’ superposition model [18]. Kriebel simulated the freak waves using the double wave superposition model consisting of a basic random wave and a linear superposition of transient wave [19]. Liu proposed a new efficient method by modifying the Longuet-Higgins model, which greatly improved the similarity between the simulated freak waves and the target spectrum structures [20]. In the early practical application, based on the numerical simulation of one- and two-dimensional space freak waves, we have calculated and analyzed the electromagnetic scattering coefficient of freak waves to study the formation mechanism, remote sensing recognition, and other related characteristics [21]. Based on the Two-Scale Method (TSM) and the Harger distribution surface, Franceschetti proposed two kinds of sea surface SAR simulator models [22]; however, this research method has great limitations because of its failure to fully consider the strong non-Gaussian statistical characteristics and velocity bunching effect of the height distribution of freak waves [23, 24]. In recent years, freak waves have attracted much attention because of their potential for serious damage to shipping and offshore structures. In practice, people pay more attention to the possibility of predicting the occurrence of freak waves, while most of the previous research on freak waves has focused on the mechanism of freak waves. The modeling flowchart that explains the main goal, the applied methods, the intermediate steps, and the obtained results is shown in Figure 1.

In this work, an improved phase modulation method for simulating freak waves is developed based on the Longuet-Higgins model and the comparative results between JONSWAP [25] spectrum and Elfouhaily [26] spectrum, and an electromagnetic scattering model of the simulated freak waves is established. The random phase correction method and its special application in freak wave simulation are presented in Section 2. According to the research conclusions, the backscattering model of the simulated freak waves is developed, and the results of the simulation were compared with the freak wave measured in Tokai of Japan in Section 3 from the perspective of B-F instability and spectral peakedness. In Section 4, the scattering calculation results are analyzed and the experimental conclusions are summarized. Moreover, a feature identification method of freak wave from its background sea surface is proposed. Finally, based on the normalization method of Z-score, the influence of the deviation coefficient on the height of the freak waves is measured in this work.

2. Numerical Simulation Model of Freak Waves Based on the Longuet-Higgins Model

2.1. Research on Random Phase Modulation Method Based on the Longuet-Higgins Model

Marine shipping and ocean engineering structures are terribly threatened by the freak waves, which leads to many maritime accidents. Because measuring the data of freak waves is difficult and the method of laboratory simulation is expensive, it is necessary to study the occurrence and evolution characteristics of freak waves through numerical simulation. Based on the Longuet-Higgins model, the simulation of normal random waves can be realized without freak waves [27]. However, the numerical calculation of the freak wave and its background waves is simulated by correcting the random phase in this work. The fixed point of the wave equation can be expressed by superposition of a large number of random cosine waves [28]:

In (1), is the time course of wave, is the sum of wave numbers, and is the distance from the simulated wave. , , , and are the amplitude, wave number of the i wave, angular frequency, and random initial phase of the composition wave. When we simulate the conventional random wave, the initial phase of the wave components is evenly distributed in . In order to simulate the freak waves in the random wave series, we need to focus the energy of the waveform. Normally, the method can be realized by adjusting the initial phase of the part composition wave. If the modulation process is not reasonable, the statistical characteristics of the numerical simulation of random wave sequence are not in conformity with the statistical characteristics of the natural wave, and the structure of the wave spectrum can be changed. Therefore, in this work, we use the following method to realize the simulation of random wave sequences of freak waves. Supposing that the freak waves are generated at the position of and at the time , we modulate to make a positive value, so when we simulate wave superposition, the wave height increases. We rewrite (1) into the following synthetic waves containing freak waves and normal waves.

We assume that the second part of the wave, , will produce the freak waves at a predetermined position; at this time, we should modulate to make a positive value. When , we make and , and we modulate , after which we can get the conclusions that and . At this time, and ; we modulate the value of as follows in this paper:When the value of is in the range of , the range of random values is When the value of is in the range of , the range of random values is When the value of is in the range of , the range of random values is When the value of is in the range of , the range of random values is

By the same way, we can deduce the value of when it satisfies the condition of .

Based on the abovementioned random phase modulation method, the research goal can be achieved. The original wave height is artificially divided into two parts: random superposition and positive superposition, and then the simulation of the freak wave is realized. The advantage of this method is that different wave height simulations can be achieved by controlling the ratio of the two parts of the superimposed wave during the experiments.

2.2. Numerical Simulation of Freak Waves and Result Analysis

The actual sea surface is often an unsteady sea surface caused by complex environments such as swells. The classic JONSWAP spectrum [12] corrects the gravity wave area on the basis of the traditional PM spectrum, so that it contains the unsteady sea spectrum. Its power spectrum is shown as follows:

In (3), is the directionless curvature spectrum, is the wavelength of the gravitational wave, is close to the PM shape spectrum parameter, and .is the peak enhancement factor, . In addition, , , , , and are wind zones in units of m. When the corresponding values of are 0.84, 1.0, and 2.0, they represent fully developed sea surface, mature sea surface, and developing sea surface, respectively.

Elfouhaily proposed a joint spectrum function based on PM spectrum, JONSWAP spectrum, and Apel spectrum, and its power spectrum is defined as follows [12]:

In (4), and , respectively, represent the low-frequency nondirectional curvature spectrum corresponding to the gravity wave and the high-frequency nondirectional curvature spectrum corresponding to the capillary wave. The low-frequency curvature spectrum satisfies the following form:


In the equations above, is the generalized P-K equilibrium zone parameter in the low-frequency wave number range, is the phase velocity corresponding to the peak of the spectrum, is the wind speed at 10 m above the sea, and is the inverse wave age. is the phase velocity, is the wave number distributed in the peak of the spectral domain, , wherein is the density of seawater, is the surface tension of seawater, and is the acceleration of gravity. Subsequently, Elfouhaily further proposed that , and based on the dimensionless parameters and , we can further determine the long-wave edge effect function .

In the comparative experiment, a simplified scattering model is used to analyze the power characteristics, , where and represent the sea surface slope in different directions, P is the sea surface slope density function, k is the wave number, R is the Fresnel scattering coefficient, and q represents the scattering vector. Based on the control variable method, the scattered power of the JONSWAP spectrum and that of the Elfouhaily spectrum are compared to normalize the power waveform peak value and the slope change rate under different wind speeds and wind conditions, wherein the receiver height is 4.5 km, the satellite elevation angle is 30°, the wind direction is 0°, the wind speed varies from 6 m/s to 20 m/s, and the fetches varies from 10 km to 19 km. The numerical results are shown in Figures 2 and 3.

Comparing the numerical results in Figure 2, it can be seen that as the wind speed becomes larger, the peaks of the JONSWAP spectrum and the Elfouhaily spectrum gradually decrease, the delay slope of the retardation gradually increases, and the effect of medium and low wind speeds is obvious. At the same time, the numerical results in Figure 3 show that, with the increase of the fetch, the peak and delay slope of the JONSWAP spectrum show regular changes, while the Elfouhaily spectrum is not sensitive to the fetch.

In this paper, the JONSWAP spectrum is used as the target spectrum [12], which means that the parameter in (1) always complies with the JONSWAP spectrum. The time series of freak waves can be simulated when the distance from the simulated wave is in Figures 4 and 5. Similarly, the space series of freak waves can be simulated when in Figure 6.

When the depth of water is , the effective wave height is , the spectral peak period is , the spectral elevation factor is , the wave number is , the modulation wave number is , the spectrum range changes from to , and the time is or ; the simulation of the time series of freak waves is illustrated in Figures 4 and 5.

According to the definition of the freak waves, the height of freak wave should meet the following conditions: , , , and , wherein is the crest height of freak waves corresponding to the horizontal line, is the effective wave height, and and are the wave heights of adjacent waves before and after the deformed wave. , , , and are characteristic parameters of freak waves [28, 29]. The characteristic statistics of wave duration are carried out using the method of positive and reverse zero-crossing counting, and the characteristic parameters of the extreme waves are shown in Tables 1 and 2.

The effective wave height is 3.76 m. Compared with the input parameter, the relative error is less than 5%. From Table 1 and Table 2, we can see that all of the parameters above meet the definition of the freak wave, and the freak wave is generated at the scheduled time, which proves the validity of this model. Wave time history spectrum and the target spectrum are compared in Figure 7, and the comparison results between wave height distribution and Rayleigh distribution are shown in Figure 8.

The results shown in Figure 6 indicate that the simulated wave spectrum keeps the structure of the target spectrum, and the spectral peak frequency is very similar to that of the target spectrum. Figure 7 shows that the normalized cumulative probability distribution of wave height of the simulated data agrees well with the Rayleigh distribution. The results show that the simulation results meet the requirements of random waveforms, and the simulation method proposed in this paper is effective.

3. Electromagnetic Scattering Calculation Model of Simulated Freak Waves

3.1. Research on Backscattering Model Based on KA Method and TSM

The Two-Scale Method (TSM) is developed on the basis of Kirchhoff approximation (KA) method adapted to the large-scale sea surface and the small perturbation method adapted to the small-scale sea surface [30, 31]. Firstly, the scattering coefficients of small-scale sea surface are calculated by the perturbation theory; secondly, the scattering coefficients of the mean sea surface are calculated considering the slope distribution of the large scale; and finally, the theory of the Two-Scale Method is used in this paper. When the incident plane is located in the space, the backscattering coefficient is calculated as follows [32]:

Here, means horizontal polarization and means vertical polarization. is the incident angle. and are the slopes of the rough surface along the x and y directions. is the probability density function satisfying the slopes of the large-scale surfaces in different directions. , , , and are the backscattering coefficient results under different polarization states. and are the backward scattering coefficients of small-scale capillary waves in the horizontal and vertical polarization. The expressions are shown below [33].

Here, HH and VV represent different polarization modes. is the direction function. is the incident wavelength, is the wave number of the incident wave, and is the incident angle. is the two-dimensional sea spectrum. The simulation of the 2-D freak wave is realized by adding the direction function on the basis of the 1D power spectrum to express the anisotropy of the energy distribution. and are the polarization amplitudes under horizontal and vertical polarization. In (8), and represent the large-scale cutoff wave number and the small-scale cutoff wave number, respectively. Among them, determines the correction of the slope probability density function in the TSM. For a given wavelength, the rough sea surface meets the condition , which constitutes the large-scale part of the Two-scale Model, whose scattering coefficient is calculated by the SPM method, and it satisfies the condition of the first order perturbation approximation. On the other hand, the small-scale rough sea surface should also contain the spatial wave number which satisfies the Bragg scattering condition of . We use the condition of to determine the large-scale cutoff wave number. However, it is necessary to satisfy the condition of and during the calculation of the scattering coefficient of small-scale rough sea surface, and the boundary threshold is directly affected by the wave number . The root mean square formula of the small-scale part is , from which we can calculate directly [34]. Corresponding to the spatial sequence simulation of freak waves based on (2), we set the parameter to 0, which represents the wave direction spectrum in .

3.2. Analysis of Electromagnetic Scattering Results

Based on the simulation conditions of freak waves mentioned above, in this work, we assumed that during the experiment the one-dimensional freak wave appears at and the adjusted deformity ratio is . The simulation results of the background wave space sequence and the freak wave space sequence are shown in Figure 6; BW means background wave, and FW means freak wave.

During the calculation of the electromagnetic scattering coefficient of the sea surface, the wind speed is 14 m/s based on the actual data, the radar operating frequency is , the incident angle is , the relative azimuth angle is , the polarization mode is VV, the sea water dielectric constant is 81, and the fetch is . Based on the parameters mentioned above, the electromagnetic scattering coefficients of the freak waves in Figure 6 are calculated, and the results are shown in Figures 9 and 10.

In Figures 9 and 10, the electromagnetic scattering characteristics of freak waves and the background waves are compared under the condition of VV polarization, the NRCS of the freak waves varies periodically with the distance of axis x, and the scattering characteristics are similar to each other. When the NRCS of background waves reaches the maximum at the crest or reaches the minimum at the trough, the NRCS of freak waves is smaller than that of background waves. The maximum difference value −47.12 dB appears at the position 100 m; the freak wave is shown in Figure 6. Comparing Figures 6 and 10, we can conclude that when the extreme wave appears in the one-dimensional simulated sea surface of background waves and freak waves, the electromagnetic scattering coefficient presents the nonsmooth transition and instability obviously. This is because when we use the Two-Scale Method to calculate the backscattering coefficient of sea surface, the sea surface is divided into large-scale gravity waves and small-scale capillary waves artificially. However, in fact, the transition of the sea surface is smooth, which indicates that the Two-Scale Method does not fully meet the physical significance. In addition, it is found that, in the vicinity of the extreme wave, the electromagnetic scattering coefficient of the background waves and the freak waves changed nonsmoothly in Figures 6 and 10, because the incident angle is in the calculation of electromagnetic scattering coefficient, which fails to fully consider shielding effect. What is most important is that the NRCS of freak waves is much smaller than that of the background waves at the position of . BW means background wave, and FW means freak wave.

According to the above simulation results, the characteristic parameters of freak waves are analyzed when the wind speed changes from 6 m/s to 20 m/s, and the experimental results are shown in Table 3 and 4. From the tables, we can find that when the wind speed is 14 m/s, the characteristic parameters of , , , and began to meet the criteria of the freak waves. At this time, the difference in NRCS between FW and BW is −12 dB; as the wind speed continues to increase, the characteristic parameters meet the criteria of freak waves stability, and the difference in NRCS also decreases steadily, which means that the difference in NRCS can also be used as a criterion for the identification of freak waves, and the decision threshold is −12 db.

3.3. Comparative Study of the Numerical Calculation Model and the Measured Data

Finally, based on the Z-score normalization method, the effect of the deviation of the deformity coefficient on the height of the freak wave is measured in this work. The formula is expressed as follows: , wherein is the standard value of the characteristic parameter of the freak wave, is the measurement level of different measurement units when calculating the average error, is the average value of different measurement types, and is the average absolute deviation. , , and the values of i are 1, 2, 3, and 4.

From Figure 11, it can be seen that the normalized malformed parameters have the same influence on the backscattering coefficient, and has the greatest influence on the NRCS of the simulated freak waves. The influences of , , and are small, and the fitting accuracy of each of them is better. Moreover, when the normalized malformed parameters are distributed in the range of , the backscattering coefficient no longer increases. At this time, the NRCS tends to be stable and reaches a critical value. At this time, and . By analyzing the experimental result in Figure 12, it can be seen that the simulated freak wave is basically consistent with the target spectrum in the distribution of scattering coefficients, and both of them have a critical smoothing phenomenon in the red region in the range of 0.52–0.67.

Based on the above numerical conclusions, the effective index of B-F instability (BFI) [34] is used in this work to test the probability of occurrence of freak waves. This index is directly related to the quasi-four-wave resonance interaction and wave surface displacement deviation from the normal distribution. The expression of the above factors is shown below [34].

Here, is the wave number of the peak frequency of the spectrum, is the zero-order distance of the wave spectrum, is a physical quantity describing the width of the wave spectrum, is the energy spectrum of the simulated freak wave surfaces, and reflects the energy distribution along different directions [35]. On June 23, 2008, Suwa Maru fishing boat carrying 20 fishermen sank in Tokai of Japan [36]. According to reports, the sea conditions were moderate and the wave height was about 2–3 meters at that time. The result of the investigation is that the fishing boat is most likely to have encountered a freak wave. A large spectral peakedness corresponds to a small spectral width. This paper uses [2] to calculate and analyze the spectral sharpness and BFI evolution over time based on the output two-dimensional ocean wave spectrum. Combined with the comparative data of the freak wave measured in Tokai of Japan (, ), the effectiveness of the numerical simulation model and scattering calculation model proposed in the work is analyzed. Among them, the wind speed at the accident site is about 11 m/s, and the wave height is about 3.50 m [37]. Figures 13 and 14 show the results of spectral peakedness and B-F instability (BFI) values of the freak wave at the point of time, respectively.

From Figure 13, it could be obviously found that the spectral peakedness is near a minimum value at the time of the accident, which is about 2.19, and the wave spectrum is wide at this time. This state indicates that the possibility of freak waves due to B-F instability at the accident site is extremely small. As shown in Figure 14, the BFI value is near the minimum value at the time of the accident, about 0.23, which is far from the condition where B-F instability easily occurs. In addition, from the perspective of time changes, the wave height is in the process of rapid changes from large to small. The measured wind speed is relatively stable in the measured data, and the wind speed has been maintained at about 11 m/s for a long time at the time of the accident, corresponding to a wave height of about 3.50 m. The above conditions indicate that the wave was not in a pure wind wave state at the time of the accident. Compared with the smooth transition coefficient of the target spectral scattering model, the freak wave simulation method based on Longuet-Higgins model is more in line with the variation law of sea surface, which proves the effectiveness of the proposed method. In addition, it is found that the average absolute deviation, used to measure the influence of the deformed parameters and backscattering coefficient on the freak waves, is more robust than the traditional method in the calculation process of NRCS. This method is more reflective on the characteristics of the data obtained by different methods. Especially in the process of normalizing outliers, the traditional deformation coefficient is very sensitive, but the Z-score normalization can achieve better robustness through the substitution method mentioned above. However, its computational complexity is higher.

4. Conclusion and Future Work

In this work, the practical problems of weakly nonlinear Longuet-Higgins model in the simulation of strongly nonlinear freak waves are fully considered. Based on the Longuet-Higgins model and the comparative experiments of wave spectrum, a modified phase modulation method for simulating freak waves is developed. The surface elevations of some wave components at the preassigned place and time are positive by modulating the corresponding random initial phases, which enhances the total surface elevation and causes a freak wave to be generated. Comparative experiments with the measured data from Tokai of Japan show that the method can make the freak waves not only occur at specified time and place but also conform to the statistical characteristics of the wave sequence and keep consistent with the target spectrum. The advantage is that the simulation of different wave heights can be achieved by controlling the ratio of superimposed waves during the experiment. In addition, this method can also make the initial phase randomly distributed between 0 and . An electromagnetic backscattering model is established based on KA method and TSM in this work, and the numerical results show that the NRCS of freak waves is much smaller than that of the background waves. Therefore, the NRCS can be used as the feature identification of freak waves, especially considering that it is difficult to obtain the characteristic parameters of , , , and in practical applications, but it is relatively easy to obtain the NRCS difference between FW and BW from SAR images. Based on the results of model data analysis, we can draw the conclusion that when the NRCS difference of SAR image is less than or equal to -12 dB, we can determine the occurrence of the freak waves. These research ideas provide a reference standard for early warning identification of freak waves in practical engineering applications.

Data Availability

No data were used to support this study.

Conflicts of Interest

The authors declare that they have no conflicts of interest regarding the publication of this paper.


This study was supported by the Shandong Key R&D Program (Soft Science) Project under Grant no. 2019RKB01058; the Shandong Provincial Natural Science Foundation, Grant no. ZR2018BD029; the Scientific Research Foundation of Shandong University of Science and Technology for Recruited Talents, Grant no. 2017RCJJ046; the SDUST Research Fund; the Training Program of the Major Research Plan of the National Natural Science Foundation of China: Research on Planting Structure Optimization and Production Management Decision Technology Based on Agricultural Big Data, Grant no. 91746104; the National Key Research and Development Program of China: Mine Internet of Things of Cloud Interaction Technology and Service Platform, Grant no. 2017YFC0804406.