Electromagnetic Scattering Analysis of the Sea Surface with Single Breaking Waves

A new electromagnetic (EM) scattering model of the sea surface with single breaking waves is proposed based on the highfrequency method in this paper. At first, realistic breaking wave sequences are obtained by solving the fluid equations which are simplified. ,en, the rough sea surface is established using the linear filtering method. A new wave model is obtained by combining breaking waves with rough sea surface using a 3D coordinate transformation. Finally, the EM scattering features of the sea surface with breaking waves are studied by using shooting and bouncing rays and the physical theory of diffraction (SBRPTD). It is found that the structure that is similar to a dihedral corner reflector between the breaking wave and rough sea surface exhibits multiple scattering, which leads to the sea-spike phenomenon that the scattering result of horizontal (HH) polarization is larger than that of vertical (VV) polarization, especially at low-grazing-angle (LGA) incidents with upwind. ,e sea-spike phenomenon is also closely related to the location of strong scattering.


Introduction
Target recognition in a complex ocean background is a focus problem that received large attention. Sea clutter influences the detection of targets in extreme sea conditions, which in turn increases the probability of false alarms. erefore, it is critical to construct a realistic breaking wave model and study its electromagnetic (EM) scattering characteristics.
LONGTANK profiles [1] laid the foundation for research on breaking waves. Waves at different times are spliced together in two-dimensional LONGTANK profiles to generate three-dimensional breaking waves [2,3]. In particular, three-dimensional breaking waves have been generated by the finite element method combined with the ocean dynamics [4]. One-dimensional time-varying breaking waves have been discussed [5]. A limited size dihedral corner reflector split structure has been used to approximate the overflow breakage [6]. Breaking waves have been determined by solving the governing equation of fluid mechanics [7,8]. However, most of the previous methods are not realistic enough and often inefficient. Moreover, they may be only employed to study small areas of sea surface containing breaking waves.
Concerning EM scattering, scholars have studied EM scattering characteristics of breaking waves based on the LONGTANK model. Holliday et al. used the method of moments (MoM) with backward and forward iteration to calculate the scattering of breaking waves with different geometrical shapes [9,10]. West and Zhiqin Zhao studied EM scattering of breaking water waves with rough faces [11]. Zhiqin Zhao and West [12] used the multilayer fast multipole method to investigate scattering characteristics of three-dimensional breaking waves. Superevents where scattering due to horizontal (HH) is larger than that due to vertical (VV) have been found in the above numerical simulations.
In this work, realistic breaking wave sequences are obtained by solving the fluid equations which are simplified. e model includes the detailed characteristics of the wave at different times and is therefore inherently efficient. Moreover, a new sea surface with the single breaking wave model is obtained by combining breaking waves with rough sea surface using a 3D coordinate transformation. It is necessary to take into account the multiple scattering produced by breaking wave itself, and the coupling effect between breaking wave and sea surface. ere are two methods for doing this. e first is the numerical method where it is suitable for calculating small objects. e second one is a high-frequrency approximation method which suites large objects. e given sea surface with a breaking wave is a large target; therefore, the shooting and bouncing rays and physical theory of diffraction (SBR-PTD) method is used here. Indeed, the SBR method effectively combines the advantages of geometrical (GO) and physical optics (PO). In particular, it is necessary to take into account the multiple scattering produced by breaking wave itself, and the coupling effect between breaking wave and sea surface.

Mathematical Model and Theoretical Formulations
2.1. ree-Dimensional Simulation of Breaking Waves. Many popular methods for modeling water surfaces work well in producing static images, but these methods only attempt to depict a surface with roughness, and they are not suitable for the realistic dynamics of the surface over time [13,14]. Moreover, traditional modeling methods based on the shallow water approximation are inefficient. erefore, an improved breaking wave model is established here building on reference [15].
In practice, breakers are normally associated with shorelines but can occur anywhere in the ocean. e major disturbing force in the open ocean is wind. In this paper, the plunging breakers are studied. ey occur when the wave encounters an abrupt transition from deep to shallow water. e base of the wave decelerates rapidly, while the top of the wave continues moving at a higher speed. With this largespeed differential, the top of the wave pitches out in front, forming a curl or tube. Moreover, we just make an assumption and establish a geometric model of the breaking waves, which is combined with the rough sea surface. e model is applied to the study of EM scattering. e process is divided into three steps: (1) e realistic breaking wave sequences are obtained by solving the fluid equations which are simplified. (2) e sea surface is modeled by the Monte Carlo method based on the Elfouhaily [16] sea spectrum. (3) A new sea surface with the single breaking wave model is obtained by combining breaking waves with the rough sea surface using a 3D coordinate transformation.
Step 1 and Step 3 will be described in detail in the following.
We begin with a vastly simplified set of equations that has been widely used for shallow water. e simplification arises from three approximations: (a) We assume that the water surface is a height field. Figure 1 shows the discrete representation of the height field in two dimensions. is, of course, has some obvious limitations. e water cannot splash and waves cannot break. However, so long as the forces on the water are sufficiently gently, the height-field assumption will not introduce error. (b) e vertical component of the velocity of the water particles can be ignored. Once again, the limitations of this assumption are fairly clear. If a disturbance creates very steep waves on the water surface, the model will cease to be accurate. (c) e horizontal component of the velocity of the water in a vertical column is assumed constant.
If there is turbulent flow or unusually high friction on the bottom, this assumption will break down.
ose assumptions lead to some obvious limitations, but the experience of hydrodynamicists suggests that they provide a rather useful model to describe breaking sea wave [15].
For simplicity, we begin with a height-field curve in two dimensions. Later, the same techniques will be extended to a height-field surface in three dimensions. Let z � h(x) be height of the water surface and let z � b(x) be the height of the ground. Let is the water depth and u(x) is the horizontal velocity of a vertical column of water. Waves where d > L/2 are called deep water and where d < L/2 are called shallow water waves [17]. In this, d is the water depth. L is the wavelength. e shallow water equations following from the above assumptions [18,19] may be written as follows: where g is the gravitational acceleration. Equation (1) expresses Newton's law F � ma, while equation (2) accounts for volume conservation. Notice that even with the above three simplifying assumptions, the resulting differential equations are still nonlinear. A further simplification which is often used is to ignore the second term in equation (1) and linearize the expression around a constant value of h. e resulting equations are then If we differentiate equation (3) with respect to y, then differentiate equation (4) with respect to t, and finally substitute for the cross-derivatives, we end up with 2 International Journal of Antennas and Propagation which is the one-dimensional wave equation with wave velocity �� � gd. In order to solve equation (5), we need to construct a discrete representation of the continuous partial-differential equation. Here, the finite-difference technique works particularly well because of the simple height-field representation. After experimenting with a number of finitedifference approximations to equations (3) and (4), the most stable version we have found is where Δy is the separation of the samples along the y direction. Putting the above two equations together, we get which provides the sought discrete approximation to equation (5). For simplicity, we use a first-order implicit method to solve ordinary differential equation (8). Let h(n) to denote h at the nth iteration and let dots denote differentiation with time. en, the first-order implicit equations may be written as Notice that the right-hand sides of these equations are evaluated at time n which corresponds to the end of the iteration rather than time n − 1 which corresponds to the beginning of the iteration. is is what makes the iteration implicit and stable. Rearranging the above expression, we arrive at By substituting equation (8) into equation (10), we obtain After these manipulations, equation (11) is still a nonlinear equation. In order to find its solution, it needs to be linearized, i.e., d is regarded as a constant, and the wave velocity is fixed as a function of y. In this linear regime, the next value of h can be calculated from previous values with the symmetric tridiagonal linear system where the matrix A is given by and the elements of A are as follows: Now, let us change the equation to where τ is phenomenologically introduced to describe damping phenomena. If τ � 0, then equation (15) reduces to equation (12), whereas if τ is between zero and one, it will make the waves damp out over time.
e effect is that observed in a viscous fluid. ere is one further subtlety of importance in the two-dimensional case. Even though equation (15) was derived from equation (6) which specifies conservation of volume, there is no guarantee that the results of the iteration will precisely conserve volume. e primary cause of departures from volume-conserving behavior is that the iteration may leave h i < b i for some index i. To compensate for this negative volume, the iteration will create excess positive volume elsewhere. While the effect is small, it can accumulate over time and create substantial drift. If the entire surface acquires a small net upwards velocity, it will very quickly create noticeable amounts of water. To combat this effect, the following simple projection appears to be adequate. After each iteration, we find the connected pieces of the fluid, and this can be done by scanning the h and b vectors in order and testing whether h i < b i . For each connected piece of the fluid, we calculate the old volume and the new volume. If the new volume is different, we distribute the difference uniformly over the samples in the connected region.
e above treatment is valid to a height field in two dimensions.
e three-dimensional case can be approximated by a series of two-dimensional equations. e basic wave equation for water in three dimensions is the same as the two-dimensional case except that the second derivative of h with respect to y is replaced by the Laplacian.
In order to solve the equations in three dimensions, we rely on the alternating-direction method [15]. e basic idea of the method is to take equation (16) and split the righthand side into the sum of two terms, one of which is independent of y, and the other is independent of x. We then divide the iteration into two subiterations. In the first, we replace the right-hand side of equation (17) with the first term, and in the second one, we replace the right-hand side of equation (17) with the second term. More specifically, in the first subiteration, we solve the equation and in the second subiteration, we solve the equation For the first subiteration, we compute the update as before on each row of the height field. For the second subiteration, we do the same for each column in the height field.
en, the breaking waves generated are combined with the rough sea surface by 3D coordinate transformation. is method will be described in detail in the following.
Rotation matrices Rx, Ry, and Rz about x, y, and z axes are defined as where ϕ, θ, and φ are angles with x, y, and z axes, respectively.
Scaling matrices are defined as where, xs, ys, and zs are scaling factors. World space matrix is defined as where xt, yt, and zt are coordinates of the breaking waves.
Let us now define (3), L is length of the breaking waves, and width is width of the breaking waves. According to [20,21], the relationship between wave height H(m) and wind speed U (m/s) is where α � U 2/3 − 12.6549 and β � U 2/3 − 6.4463. Wave steepness S [22] is defined as the ratio of H/L, i.e., where S is the wave steepness. We can get the wavelength L (m) according to (25) e smoothing factor is defined as where amplitude � 1. If nPtLen < 0, then nPtLen � 0, whereas if nPtLen > 1e 9 , then nPtLen � 1e 9 . Convert positional coordinates from local to global coordinates, i.e., where bwc are the coordinates of the sea surface with single breaking waves. e sea surface with single breaking waves is shown in Figure 2. e evolution process describes the breaking waves from the initial moment to gradual curl. It is worth mentioning that the wave is a model that changes all the time, and Figure 2 is only a realistic description of the state of the wave probably. e time interval is taken as Δt � 0.05 s. e size of rough sea surface is 5 m × 5 m with wind speed U � 5 m/s. S � 1/3, and the width of the breaking waves is 3 m. e height and length of breaking waves can be calculated by equations (25) and (27), respectively. erefore, the location of the breaking wave can be determined. e aim of this paper is to study EM scattering characteristics in the presence of the sea surface with single breaking waves, so this paper only studies the breaking waves before involution.

SBR-PTD Method.
e SBR method is a hybrid method of GO combined with PO, and the PTD method can calculate edge diffraction. Calculation accuracy can be improved using PTD together with SBR. e SBR-PTD method uses GO to determine the triangles illuminated by the incident wave and each order of reflected wave. Moreover, recording the IDs of each triangle's three neighbors is also useful in the calculation of PTD fields. e tracing process may be understood by looking at the case depicted in Figure 3, where a ray from triangle patch 0 is shot toward a target, which consists of 16 triangle patches. e following steps should be carried out to determine which triangles are illuminated by the ray from triangle patch 0. At first, forward ray tracing is performed to detect one of these illuminated triangles (i.e., triangle 7). en, after finding illuminated triangle 7, it is necessary to determine whether the adjacent triangles (i.e., triangles 6, 8, and 13) are illuminated as well (e.g., triangle 13). Backward ray tracing is used to determine whether triangle 13 intersects triangle patch 0 and is not occluded. If so, triangle 13 is illuminated; otherwise, it is discarded. en, backward ray tracing is used to determine whether triangle 6 and triangle 8 are illuminated or not, and the procedure continues until some triangles outside the illuminated area occur. Moreover, a binary tree is built to accelerate the process of ray tracing.
For each triangle patch, the scattered far field is obtained by vector-summing each order of the PO fields from the triangle surface and the PTD fields from the three edges of this triangle [23].
e electric field, approximated by PO integral, may be the written as follows: where k � ω(μ 0 ε 0 ) 1/2 , η � (ε 0 /μ 0 ) 1/2 , s is the unit vector aligned along the scattering direction, R is the distance from the specular point r → ' to the observation point, and J → is the current density on the triangle surface. e diffraction field from each edge of a triangle patch is approximated by the PTD method: International Journal of Antennas and Propagation 5 where I e and I m are as follows: where E → inc and H → inc denote the incident electric field strength of patches' edges; β i/s stands for the angle between the edge and incident or scattered ray; and D PTD e , D PTD m , and D PTD em represent the PTD diffraction coefficients [24], which are closely related to the angle between two adjacent triangles whose IDs have already been stored in the computer memory.

Verification of the Effectiveness of the SBR-PTD Method.
Back-scattering RCS of the sea surface with single breaking waves both in t � △t stage and t � 11△t stage is shown in Figure 4. e frequency of the incident wave is f � 3 GHz. erefore, the dielectric parameter is ε � (70.6059, 39.7343) from the Debye model [25]. e incident angles are θ i � −90°∼ 90°, φ i � 0°. e size of rough sea surface is 1 m × 1 m with wind speed U � 5 m/s. It may be readily seen that the result of SBR-PTD agrees well with the result of MoM, which verified the effectiveness of the SBR-PTD method. Note that θ i > 0°r epresents downwind and θ i < 0°r epresents upwind, and this is true for all of the following examples.

A Comparison of the Bistatic Scattering Coefficient with and without Breaking Waves.
A locally breaking wave resembling a dihedral corner reflector is shown in Figure 5. A comparison of the bistatic scattering coefficient with and without breaking waves is shown in Figure 6. In this, the t � 12△t stage of the sea surface with single breaking waves is selected. e frequency of the incident wave is f � 10 GHz. erefore, the dielectric parameter is ε � (55.8866, 37.6777) from the Debye model. e incident angles are θ i � 30°, φ i � 0°. e size of rough sea surface is 5 m × 5 m with wind speed U � 5 m/s. In Figure 6, for the sea surface with single breaking waves, there is an obvious enhancement in the backward direction. is is caused by the fact that the structure between the breaking waves and the rough sea surface (i.e. Figure 5) is similar to a dihedral corner reflector, and thus it exhibits multiple scattering, resulting in the change of scattering field near the backward nonspecular direction which is small, while the sea surface scattering has a dominant contribution at the specular direction.

Back-Scattering RCS of Sea Surface with Single Breaking
Waves for Different Stages. To numerically demonstrate the sea-spike phenomenon of breaking waves, t � △t, t � 5△t, International Journal of Antennas and Propagation t � 8△t, and t � 11△t, which denote four different time evolution histories, are chosen to analyze the EM scattering mechanism. Back-scattering RCS of the sea surface with single breaking waves in the t � △t, t � 5△t, t � 8△t, and t � 11△t stages is shown in Figure 7. e frequency of the incident wave is f � 15 GHz. erefore, the dielectric parameter is ε � (43.8890, 39.1149) from the Debye model. e incident angles are θ i � −90°∼ 90°, φ i � 0°. e size of rough sea surface is 5 m × 5 m with wind speed U � 5 m/s. In Figure 7, for t � △t stage of the sea surface with single breaking waves, it is at the very beginning of the plume forming. e multipath effects are not serious. HH is almost equal to VV. For the t � 5△t, t � 8△t, and t � 11△t stages of the sea surface with single breaking waves, the plume forms. HH is larger than VV at the grazing angles of less than 30°, and the sea-spike phenomenon has occurred. is is caused by the fact that the structure is similar to a dihedral corner reflector between the breaking    International Journal of Antennas and Propagation waves and the rough sea surface (i.e., Figure 5), and thus it exhibits multiple scattering which, in turn, enhances the Brewster effect, which will cause the reflection coefficient to decrease at VV polarization, and the VV multipath is greatly attenuated.

Back-Scattering Field Intensity Distribution of the Sea
Surface with Breaking Waves. e back-scattering field intensity distribution for different incidence angles of each triangle patch in the t � 12△t stage of the sea surface with single breaking waves is shown in Figure 8. e frequency of the incident wave is f � 10 GHz and the corresponding dielectric parameter is ε � (55.8866, 37.6777). e size of rough sea surface is 5 m × 5 m and the wind speed U � 5 m/s. e incident angles are 30°(downwind/upwind), 85°( downwind), and 87°(upwind), respectively. \scale 90%For θ i � 30°(downwind/upwind), most of triangle patches on the model are illuminated. Only the inside of the breaking waves is obscured. At this time, the scattering of the rough sea surface is dominant. For θ i � 85°/87°(downwind/upwind), a few triangle patches on the model are illuminated by EM waves. Most of the surface is partially obscured by the breaking waves of the model at low grazing angle (LGA). At this time, the scattering of the breaking waves is dominant. In Figure 8, for moderate incident angles, the main scattering comes from the rough sea surface of the model. For LGA, the main scattering comes from the breaking waves of the model. Moreover, a few triangles patches with warm colors have a stronger scattering effect compared to other triangle patches with cool colors, which reflects the strong scattering location of the model, i.e., the structure is similar to a dihedral corner reflector formed between the breaking waves and the rough sea surface (i.e., Figure 5).

Influence of Incidence Angle on EM Scattering of Sea
Surface with Single Breaking Waves. Back-scattering RCS for different incidence angles in the t � 12△t stage of the sea surface with single breaking waves is shown in Figure 9. e frequency of incident wave is f � 10 GHz, and the dielectric parameter ε � (55.8866, 37.6777) from the Debye model. e size of the rough sea surface is 5 m × 5 m with wind speed U � 5 m/s. Incidence angles are 25°(downwind/upwind), 45°(downwind/upwind), and 87°(downwind/upwind), respectively.
In Figure 9, for θ i � 25°(downwind/upwind), HH is larger than VV at some angles. However, overall, HH is almost equal to VV. For θ i � 45°(downwind/upwind), HH is larger than VV at some angles, and the sea-spike phenomenon occurs. For θ i � 87°(downwind/upwind), HH is larger than VV, even at some moderate incident angles, and this phenomenon is even more evident in the upwind. erefore, in Figures 8 and 9, it can be inferred that the main reason of sea-spike occurrence is due to the multiple scattering of dihedral corner reflector between the breaking waves and the rough sea surface (i.e. Figure 5), which will enhance the Brewster effect, and the VV multipath is greatly attenuated.

Scattering Echo Time Series.
e time behavior of wave back-scattering RCS for different incident angles and HH polarization is shown in Figure 10. In Figure 10, the incidence angles of downwind and upwind are set as 30°, 45°, 60°, 85°, 86°, and 87°, respectively, and 64 time points of the model are selected, which describe the breaking waves gradually curling. e incidence wave frequencies are set to f � 1.5 GHz, f � 3 GHz, and f � 10 GHz, respectively. e time interval is taken as Δt � 0.05 s. e size of the rough sea surface is 5 m × 5 m with wind speed U � 5 m/s (dielectric constants of sea wave are obtained with the Debye model).
In Figure 10, one can see that for increasing frequency, the back-scattering RCS changes with the angle (downwind/upwind) are more and more pronounced. For a certain time point, the back-scattering RCS of the model decreases gradually as angle increases in the downwind case. However, the change of the back-scattering RCS is not regular when the angle increases at the upwind incident. Moreover, for a certain angle, back-scattering RCS at each angle is gradually increasing as time goes by at the upwind incident. e main reason for this behavior is again the fact that the structure is similar to a dihedral corner reflector between the wave and the sea surface. e structure is gradually formed, and this enhances the Brewster effect.

Conclusions
In this work, a sea surface with single breaking waves which contains realistic details has been established. Moreover, the SBR-PTD method is firstly applied to study the EM scattering of breaking waves. e existence of superevents where the scattering of HH polarization is larger than that of VV polarization has been confirmed by studying the scattering field of the sea surface with the single breaking wave model. ose superevents are more likely to occur at LGA with upwind. Moreover, the strong scattering localization of the model is also analyzed. e scattering of the sea with multiple breaking waves will be studied in the future.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that they have no conflicts of interest.