Improvement of SSA Approach for Numerical Simulation of Sea Surface Scattering at High Microwave Bands

: Small slope approximation (SSA) is a widely accepted approach in sea surface electromagnetic (EM) scattering studies. Nevertheless, the spatial sample interval used for sea surface should be around or even smaller than one-eighth of the incident wavelength to ensure EM scattering calculation accuracy, which requires a huge amount of computation, creating an obstacle to scattering numerical simulation, especially for high microwave band incident waves and large sea surface scenes. In this paper, a novel realization approach for SSA is proposed to signiﬁcantly decrease the computation demands and computer memory requirements in sea surface scattering simulation. First, the sea surface is decomposed into two scales, and each scale has its own spatial sample interval. Then, the inclination state of the large-scale sea surface is determined under a speciﬁc wind speed. After that, scattering calculations of a typical surface cell with a ﬁnely sampled structure are completed and saved in all possible situations. Finally, scattering results for all the cells of a concrete sea surface are extracted from the saved cell scattering data base. From the different kinds of scattering result comparisons, it is demonstrated that this novel SSA realization approach can attain almost similar scattering results to exact SSA. This approach can be broadly applied in composite scattering studies, and remote sense imaging simulation of large sea surfaces with multiple targets.


Introduction
The ocean is playing an increasingly important role in politics, economy, culture, and other areas worldwide. Accordingly, the studies on electromagnetic (EM) scattering simulations of marine environments are of considerable value, with applications in marine remote sensing, as well as civil and military fields [1][2][3][4][5][6]. Due to the constraints imposed by physical and financial resources, experimental data are limited in terms of the variation in sea conditions, radar parameters and viewing geometry, which has restricted the development of modern radar techniques, such as high-resolution imaging and target detections. As such, sea scattering simulation based on the EM scattering model has become a viable alternative, due to its low cost and easy implementation.
However, EM scattering faces challenges when applied to rough surfaces, so a number of approaches have been introduced in previous studies to address this. These approaches can be divided into two categories: numerical methods and theoretically approximate approaches. Numerical methods, such as method of moment (MoM) [7] and finite-difference time-domain method (FDTD) [8], usually provide the most accurate results. However, the acquisition of these accurate results comes at the cost of considerable computation requirements, which constrain the ability of numerical methods to be applied to electrically small objects or one-dimensional problems. As a comparison, theoretically approximate approaches [9,10], like the Kirchhoff approximation (KA) method, small perturbation method (SPM), two-scale model (TSM), and small slope approximation (SSA) [11], usually accomplish the calculation with their corresponding assumptions, with which the original problems can be greatly simplified, resulting in a considerable decrease in computation demand. Armed with these approaches, EM scattering of two-dimensional (2D) surfaces and some other practical problems can be solved efficiently with precise or acceptable results. Comparisons of the effectiveness and validity of these models have been discussed in detail [12,13]. Among the theoretically approximate approaches, SSA is an effective candidate that bridges the gap between the KA and SPM models, and can be degenerated when appropriate. Moreover, as SSA avoids the empirically scale partition in the TSM model, is more accurate, and has larger application scope, it has been widely applied to the practical rough surface scattering problems, and especially in sea surface scattering studies, where SSA has been proven to be reliable compared with the other methods and experiment data [14][15][16][17]. Additionally, the SSA method has been developed into a deterministic model, which is favorable for some other application areas, such as composite scattering problems for rough surfaces with targets [18], synthetic aperture radar (SAR) image simulation of maritime scenes, and dynamic scattering of the sea surface [19,20].
However, for SSA in the EM scattering calculation, the spatial sample interval for the sea surface should be around or even smaller than one-eighth of the incident wavelength to ensure accuracy of the integral operation in the SSA formula. On the other hand, the size of the simulated sea surface must be larger than its dominant wavelength to reflect the time-variation and modulation effects of the maritime scenes. This directly results in the huge number of sample points and the tremendous amount of computation required for corresponding simulation applications. Particularly, for the X-band or higher microwave bands, the spatial interval can reach the millimeter level, and the total number of sample points will be in the billions or even more [20]. To solve this problem, some assumptions or techniques have been proposed to simplify and implement the numerical simulation. Li et al [21] replaced the capillary wave with a specific sinusoidal wave according to the Bragg scattering mechanism, to simplify the SSA. Jiang et al. [22,23] applied the spectrum decomposition method to split the sea profile into two kinds of scales and filled the large scales with one small scale sample, which provided a solution for the storage issues encountered with high microwave band and large surface area. Wang and Xu [24] decomposed the spectrum of a 2D sea surface into multiple blocks and combined computer memory and external storage to accomplish the Doppler simulation of a 2D sea surface up to Ku-band. Regardless, the above approaches either lose the ability to accurately describe the sea surface or do not improve the computation efficiency during the numerical simulation.
In this paper, a novel realization approach for SSA is proposed to significantly decrease the computation amount and computer memory requirements for sea surface scattering simulation. The architecture of the proposed method is displayed in Figure 1. First, the sea surface is decomposed into two scales, and each scale has its own spatial sample interval. Then, the inclination state of the large-scale sea surface is determined under a specific wind speed. After that, scattering calculations for a typical surface cell with finely sampled structure are carried out and saved in all possible inclination situations. Finally, the scattering results of all the cells from a concrete sea surface are extracted from the saved cell scattering dataset, and the total field of the whole sea surface is obtained by summing the individual cell results.
The remainder of this paper is organized as follows: Section 1 briefly presents the SSA theory for sea surface scattering under the tapered incident wave; Section 2 outlines the detailed process for the novel implementation of SSA, including the sea surface decomposition, the scattering calculation of a typical surface cell in all possible inclination situations, and the interpolation and extraction operations. Section 3 provides the numerical simulation of backscattering scattering, bistatic scattering, and comparisons to display the feasibility and accuracy of the proposed approach. Finally, Section 4 concludes this paper. scattering, and comparisons to display the feasibility and accuracy of the proposed approach. Finally, Section 5 concludes this paper.

SSA under Tapered Incident Wave
In SSA, the geometrical configuration adopted to resolve the wave-scattering problem from the 2-D sea surface is illustrated in Figure 2, where we consider the sea surface Im 0 q ≥ ) are horizontal and vertical projections of the incident wave vector, respectively [9]. λ is the wavelength in the vacuum. Before the EM scattering calculations, an important problem is that the sea surface to be simulated is of limited size, which means the surface current is forced to be zero outside the surface. This abrupt change in surface current from nonzero to zero can cause an artificial reflection from the boundaries. To avoid this, the tapered incident wave is usually used so that the intensity of incident wave can gradually decay to zero in a Gaussian manner from inside to the boundaries [25,26]. The tapered incident wave can be expressed as:

SSA under Tapered Incident Wave
In SSA, the geometrical configuration adopted to resolve the wave-scattering problem from the 2-D sea surface is illustrated in Figure 2, where we consider the sea surface z = h(r), with r(x, y) = xx + yŷ, between two homogenous half-spaces with permittivity ε 1 (upper half-space, z > 0) and ε 2 (lower half-space, z < 0) [11]. The time factor of the incident wave is assumed to be exp(−iωt). Parameters θ i and θ s are incident and scattering elevation angles, respectively, and φ i and φ s are the incident and scattering azimuth angles, respectively. The ripple surface is illuminated by a monochromatic plane wave from the upper half-space, whose incident direction is determined by wave vector where the vector k 0 = k 0xx + k 0yŷ and q 0 = K 2 0 − k 2 0 (K 0 = 2π √ ε 1 µ 1 /λ, Imq 0 ≥ 0) are horizontal and vertical projections of the incident wave vector, respectively [9]. λ is the wavelength in the vacuum. scattering, and comparisons to display the feasibility and accuracy of the proposed approach. Finally, Section 5 concludes this paper.

SSA under Tapered Incident Wave
In SSA, the geometrical configuration adopted to resolve the wave-scattering problem from the 2-D sea surface is illustrated in Figure 2, where we consider the sea surface Im 0 q ≥ ) are horizontal and vertical projections of the incident wave vector, respectively [9]. λ is the wavelength in the vacuum. Before the EM scattering calculations, an important problem is that the sea surface to be simulated is of limited size, which means the surface current is forced to be zero outside the surface. This abrupt change in surface current from nonzero to zero can cause an artificial reflection from the boundaries. To avoid this, the tapered incident wave is usually used so that the intensity of incident wave can gradually decay to zero in a Gaussian manner from inside to the boundaries [25,26]. The tapered incident wave can be expressed as:  Before the EM scattering calculations, an important problem is that the sea surface to be simulated is of limited size, which means the surface current is forced to be zero outside the surface. This abrupt change in surface current from nonzero to zero can cause an artificial reflection from the boundaries. To avoid this, the tapered incident wave is usually used so that the intensity of incident wave can gradually decay to zero in a Gaussian manner from inside to the boundaries [25,26]. The tapered incident wave can be expressed as: where where R = (r, z) = (x, y, z), t = t x + t y , γ is the parameter that controls the tapering of the incident wave, andê i α i (k 0 ) is the unit vector defining the polarization of the incident wave, which is given for α i = 1, 2 by the following expressions: Analogously, the unit vector defining the polarization of the scattered wave,ê i α s (k)(α s = 1, 2), has a similar expression:ê s 1 (k) = −(k 2ẑ + qk)/Kk andê s 2 (k) = (ẑ × k)/k.
Then, we consider the scattered wave in the upper half-surface for directions determined by wave vector K = k − qẑ, with k = k xx + k yŷ and q = K 2 − k 2 (K = 2π √ ε 2 µ 2 /λ,Imq 0 ≥ 0). The scattered wave can be defined as: With the four coefficients given by α i = 1, 2 and α s = 1, 2, the scattering process is determined by the scattering amplitude matrix S(k, k 0 ): From the SSA theory [9], the first-order scattering amplitude matrix has the following form: where B 1 is a 2 × 2 matrix containing information about the scattering process for different polarizations [11].
With the scattering amplitude matrix acquired, the scattered field and normalized scattering cross section (NRCS) for one surface sample can be obtained using the stationary approximation in Equation (8): (12) where P inc is the incident wave power received by the rough surface. Moreover, for different surface profiles, the averaged NRCS for randomly surface is:

Sea Surface Decomposition and Synthesis
To perform sea surface scattering simulations with the SSA for high microwave bands and large maritime scenes, such as backscattering and bistatic scattering coefficient predictions, composite scattering calculation of the sea surface with targets, and imaging simulation in remote sensing applications, the spatial sample interval of the sea surface is usually restricted to a value around or even smaller than one-eighth of the incident wavelength, which should be broken up to improve computation efficiency and decrease the memory requirements. Accordingly, a novel approach for the SSA implementation is proposed in this study, which includes the sea surface decomposition, typical cell calculation for all possible inclination situations, extraction, and synthesis for total field derivation. The following section provided a detailed outline of the process. With regard to dynamic sea surface generation, the 2D surface profile can be obtained by the widely used linear filter method with a fast Fourier transform (FFT) operation as follows: where r = (x, y), k s = (k s x ,k s y ) ; t is the time parameter; L x and L y are the lengths of the rectangular sea surface along the x-and y-direction, respectively; and ς(k s , t) is the sea surface spectrum-related random variable, which can be expressed as: where W(k s ) is the sea surface spectrum that is mainly determined by wind speed and wind direction; γ(k s ) is the complex Gaussian random series with zero mean and unite variance; * is the sign of complex conjugate; and ω(k s ) is the dispersion factor of the sea surface wave.
The above description provides a brief overview of the sea surface generation process. However, in the proposed realization of the SSA, the surface profile to be simulated should be decomposed into the two following types: where k cut is the truncation wave number used to classify the long and short waves; L l x , L l y , L s x , L s y are the total lengths for different profile types and dimensions; and d l x , d l y , d s x , d s y are the spatial sampling intervals for different profile types and dimensions. Using this classification two kinds of profiles are generated, large surface and small surface, each with its own simulation area and spatial sampling interval. Specifically, the large surface profile corresponds to a large simulated maritime scene, but with relatively much larger sampling intervals than one-eighth of the incident wavelength. The small surface profile determines the detailed structure of the final sea surface profile and the scattering results, of which the sampling interval is fine enough and usually smaller than one-eighth of the incident wavelength. More importantly, the size of the small surface should be properly selected to balance the computation efficiency and the decorrelation of the short sea waves, because the small surface is used to replace all the detailed structures in the entire maritime scene, which will significantly reduce the computation amount and memory requirement in the following scattering calculation. Figure 3 shows the schematic diagram of sea surface synthesis. In this figure, there are three kinds of grids: picture (a) corresponds to the whole maritime scene and is meshed by cells of the relatively large size compared with wavelength; picture (c) corresponds to the large roughness profile of a specific cell and is meshed by the large surface intervals mentioned in Equation (16); and picture (d) corresponds to the small roughness profile of a specific cell but is meshed by the small surface intervals mentioned in Equation (16). In the sea surface simulation, different kinds of rough profiles are first generated, and then the whole sea surface is obtained after synthesis processing. In other words, the final sea surface profile is achieved based on the large roughness surface profile, which is entirely paved and refined by the finely sampled small roughness surface profile without overlap.

Scattering Calculation of Typical Cell for All Possible Inclination Situations
The above section generates the contour of the sea surface profile. However, this sea surface generation is only a final profile display of the sea surface to be calculated for the EM scattering simulation. The scattering calculation is actually not performed on the whole sea surface but only on a typical cell of the fine sea surface profile and the same size as the small surface. To calculate the scattering results of the whole sea surface by only using the results from the cell surface, the scattering of the meshed facets of the cell surface will be calculated under all possible inclinations. Figure 4 shows the establishment of the dataset of a typical cell for the EM scattering calculation of all possible inclinations. First, the small surface is meshed into square grids with the sampling intervals of ,

Scattering Calculation of Typical Cell for All Possible Inclination Situations
The above section generates the contour of the sea surface profile. However, this sea surface generation is only a final profile display of the sea surface to be calculated for the EM scattering simulation. The scattering calculation is actually not performed on the whole sea surface but only on a typical cell of the fine sea surface profile and the same size as the small surface. To calculate the scattering results of the whole sea surface by only using the results from the cell surface, the scattering of the meshed facets of the cell surface will be calculated under all possible inclinations. Figure 4 shows the establishment of the dataset of a typical cell for the EM scattering calculation of all possible inclinations. First, the small surface is meshed into square grids with the sampling intervals of d l x , d l y in the x-and y-dimension, respectively. Then, the square grids are divided into triangles to carry out the inclining modulation of the large-scale sea surface waves. A large sea surface is generated with the same sampling intervals of d l x , d l y and then divided into triangles as above. The variation ranges of the normal vector inclination angles, θ and φ versus the xand y-axis, respectively, are determined for all the triangles of the large sea surface. The variation ranges of θ and φ are sampled and allocated into an inclination matrix with all possible variation pairs (θ m , φ n ). Next, one facet from the small surface cell is moved to the origin of coordinates and modulated by an inclination pair (θ m , φ n ). Then, the EM scattering calculation is performed for this specific facet and inclination by SSA. When loops on the facets and inclination matrix are implemented, the scattering results for all typical cell facets and possible inclinations are stored, and the important dataset is established.

Scattering Calculation of Typical Cell for All Possible Inclination Situations
The above section generates the contour of the sea surface profile. However, this sea surface generation is only a final profile display of the sea surface to be calculated for the EM scattering simulation. The scattering calculation is actually not performed on the whole sea surface but only on a typical cell of the fine sea surface profile and the same size as the small surface. To calculate the scattering results of the whole sea surface by only using the results from the cell surface, the scattering of the meshed facets of the cell surface will be calculated under all possible inclinations.   As described above, the data set is a three-dimensional (3D) array that stores the complex scattering field of each typical cell facet for all possible inclination angles. In other words, this set can generate the complex scattering field of all the triangle facets from a particular large sea surface, except for a phase factor concerning the position of the facet. Figure 5 provides examples of the dataset results of a particular facet for a backscattering case, where the sampling interval d l x , d l y is 0.1 m, the incident wave frequency is 14 GHz, the wind speed in the sea surface simulation is 5 m/s, and the wind direction is along the x-axis. Figure 5a,b show the amplitude value distribution of the scattering field varying with the inclination angles θ and φ at the incident angle of θ i = 20 • and θ i = 60 • , respectively. As suggested by the figures, the variation in the ranges of the parameters θ and φ are 61 to 117 and 67 to 113 • , respectively, which appears different for the maximum inclination state for different directions due to the specification of the wind direction and sea surface profile. Amplitude levels and the distribution features for different incident angles are significantly different, due to the incident angle and the facet geometry. This finding is consistent with the characteristics of rough surface scattering. Notably, textures in both figures changed smoothly, which lays the foundation for the following interpolation operation and the entire proposed approach.
Remote Sens. 2018, 10, x FOR PEER REVIEW 7 of 14 As described above, the data set is a three-dimensional (3D) array that stores the complex scattering field of each typical cell facet for all possible inclination angles. In other words, this set can generate the complex scattering field of all the triangle facets from a particular large sea surface, except for a phase factor concerning the position of the facet. Figure 5

Cell Scattering Extraction and Total Field Synthesis
Once the dataset was established of the typical cell for EM scattering calculation of the all possible inclinations, scattering results of the entire large sea surface were derived through extraction and synthesis operations. However, the inclination angles of the meshed triangles from the large sea surfaces varied continuously, but the dataset built above for inclinations was a discrete 3D array and usually had a large angle sampling interval to improve the computation efficiency. To avoid the loss in computation accuracy, an interpolation manipulation was carried out before the total field derivation. Figure 6 provides an example of the comparison of the lines from the original dataset and after interpolation, which corresponds to the line from the Figure 5b at 60 θ =°. As displayed in the figure, the original dataset was calculated at the angle sampling interval of 1° for both θ and φ , whereas the dataset after interpolation was much more refined and varied with the interval of 0.1°. When the refined dataset for possible inclinations was acquired, the total field derivation of the large sea surface was accomplished using the following operations.

Cell Scattering Extraction and Total Field Synthesis
Once the dataset was established of the typical cell for EM scattering calculation of the all possible inclinations, scattering results of the entire large sea surface were derived through extraction and synthesis operations. However, the inclination angles of the meshed triangles from the large sea surfaces varied continuously, but the dataset built above for inclinations was a discrete 3D array and usually had a large angle sampling interval to improve the computation efficiency. To avoid the loss in computation accuracy, an interpolation manipulation was carried out before the total field derivation. Figure 6 provides an example of the comparison of the lines from the original dataset and after interpolation, which corresponds to the line from the Figure 5b at θ = 60 • . As displayed in the figure, the original dataset was calculated at the angle sampling interval of 1 • for both θ and φ, whereas the dataset after interpolation was much more refined and varied with the interval of 0.1 • . When the refined dataset for possible inclinations was acquired, the total field derivation of the large sea surface was accomplished using the following operations. Remote Sens. 2018, 10, x FOR PEER REVIEW 8 of 14 Figure 6. Comparison of the lines from original data set and that after interpolation. Figure 7 shows the implementation of the cell profiles scattering extraction and synthesis for the total field derivation. During this process, the large sea surface for the EM scattering calculation was first generated with the same sampling intervals of , l l x y d d and divided into triangles as was performed for the typical cell. Then, the sea surface was divided into cells the same size as the typical cell mentioned above. For each divided cell, its corresponding 3D inclination array was then calculated, which contained all the inclination information of the included triangles, namely the θ and φ values for the x-and y-axis, respectively. The continuous distribution of the values of the calculated 3D inclination array caused mismatching between the inclination array and the discrete dataset. Accordingly, nearest value judgment was carried out among the discrete intervals of 0.1°. With the new evaluated 3D inclination array, the scattering field for each cell triangle was extracted from the refined dataset. So far, the scattering result of a particular cell from the large sea surface was derived from its corresponding inclination array and the refined dataset. This procedure was repeated for all the divided cells from the large sea surface, including completing the scattering field for every triangle from the sea surface. The total field corresponding to the whole surface was derived with high efficiency. Figure 7. Implementation of the cell profile scattering extraction and synthesis for total field derivation.
From the above description, the computation efficiency of the scattering calculation of large sea surface was significantly improved due to the application of the dataset, which avoids the repetitive scattering computation of triangles of the same inclination from the large sea surface. In establishing the dataset and taking the parameters in the above section as an example, 57 × 47 inclination situations are possible, and for each situation, a 1 × 1 m 2 square area was calculated by the SSA. Accordingly, the total computation amount is equivalent to the EM scattering computation amount of a sea surface with an area of 57 × 47 m 2 . Importantly, this total computation amount is independent of the scale of the maritime scene to be simulated, and the efficiency improvement  Figure 7 shows the implementation of the cell profiles scattering extraction and synthesis for the total field derivation. During this process, the large sea surface for the EM scattering calculation was first generated with the same sampling intervals of d l x , d l y and divided into triangles as was performed for the typical cell. Then, the sea surface was divided into cells the same size as the typical cell mentioned above. For each divided cell, its corresponding 3D inclination array was then calculated, which contained all the inclination information of the included triangles, namely the θ and φ values for the xand y-axis, respectively. The continuous distribution of the values of the calculated 3D inclination array caused mismatching between the inclination array and the discrete dataset. Accordingly, nearest value judgment was carried out among the discrete intervals of 0.1 • . With the new evaluated 3D inclination array, the scattering field for each cell triangle was extracted from the refined dataset. So far, the scattering result of a particular cell from the large sea surface was derived from its corresponding inclination array and the refined dataset. This procedure was repeated for all the divided cells from the large sea surface, including completing the scattering field for every triangle from the sea surface. The total field corresponding to the whole surface was derived with high efficiency.
Remote Sens. 2018, 10, x FOR PEER REVIEW 8 of 14 Figure 6. Comparison of the lines from original data set and that after interpolation. Figure 7 shows the implementation of the cell profiles scattering extraction and synthesis for the total field derivation. During this process, the large sea surface for the EM scattering calculation was first generated with the same sampling intervals of , l l x y d d and divided into triangles as was performed for the typical cell. Then, the sea surface was divided into cells the same size as the typical cell mentioned above. For each divided cell, its corresponding 3D inclination array was then calculated, which contained all the inclination information of the included triangles, namely the θ and φ values for the x-and y-axis, respectively. The continuous distribution of the values of the calculated 3D inclination array caused mismatching between the inclination array and the discrete dataset. Accordingly, nearest value judgment was carried out among the discrete intervals of 0.1°. With the new evaluated 3D inclination array, the scattering field for each cell triangle was extracted from the refined dataset. So far, the scattering result of a particular cell from the large sea surface was derived from its corresponding inclination array and the refined dataset. This procedure was repeated for all the divided cells from the large sea surface, including completing the scattering field for every triangle from the sea surface. The total field corresponding to the whole surface was derived with high efficiency. Figure 7. Implementation of the cell profile scattering extraction and synthesis for total field derivation.
From the above description, the computation efficiency of the scattering calculation of large sea surface was significantly improved due to the application of the dataset, which avoids the repetitive scattering computation of triangles of the same inclination from the large sea surface. In establishing the dataset and taking the parameters in the above section as an example, 57 × 47 inclination situations are possible, and for each situation, a 1 × 1 m 2 square area was calculated by the SSA. Accordingly, the total computation amount is equivalent to the EM scattering computation amount of a sea surface with an area of 57 × 47 m 2 . Importantly, this total computation amount is independent of the scale of the maritime scene to be simulated, and the efficiency improvement From the above description, the computation efficiency of the scattering calculation of large sea surface was significantly improved due to the application of the dataset, which avoids the repetitive scattering computation of triangles of the same inclination from the large sea surface. In establishing the dataset and taking the parameters in the above section as an example, 57 × 47 inclination situations are possible, and for each situation, a 1 × 1 m 2 square area was calculated by the SSA. Accordingly, the total computation amount is equivalent to the EM scattering computation amount of a sea surface with an area of 57 × 47 m 2 . Importantly, this total computation amount is independent of the scale of the maritime scene to be simulated, and the efficiency improvement increases with increasing Remote Sens. 2018, 10, 1021 9 of 14 scale of the simulated sea surface. This property provides the proposed method with an outstanding advantage in large sea surface EM scattering simulation applications.

Numerical Simulations and Discussion
In this section, some numerical simulations are illustrated for sea surface scattering by applying the proposed method, including back and bi-static scattering coefficients varying with incident or scattering angles. Moreover, in order to demonstrate the validity of the proposed method it is compared to statistical SSA.
First, the backscattering simulation of sea surface was carried out at the Ku-band (14.0 GHz) with a wind speed of 5 m/s, which was chosen for the comparison with the statistical SSA results reported previously [27]. In this simulation, the large mesh grids were set to 0.1 × 0.1 m 2 , and the small grids were 0.001 × 0.001 m 2 to accurately estimate the Bragg scattering contribution. During the scattering field dataset construction, the typical cell included 100 large grid cells of the all possible inclinations for EM scattering calculation. In the implementation process, there were 1024 sampling points in both x and y directions, and the simulation area of the sea surface was 102.4 × 102.4 m 2 . The Elfouhaily model [28] was adopted for the sea roughness spectrum, which was developed based on available field and wave-tank measurements. The Elfouhaily model is supported by strong physical arguments, contrary to other spectra that are mostly empirical. The relative dielectric constant of the sea water was calculated using the Klein dielectric constant model [29] at 20 • C and 35 ‰ salinity. With a 5 m/s wind speed, the variation ranges of the normal vector inclination angles, θ and φ versus the xand y-axis, were 61-117 and 67-113 • , respectively. After the dataset of the typical cell of the all possible inclinations was established for EM scattering calculation, interpolation was performed to obtain a much more delicate variation description with the interval of 0.1 • . Finally, the scattering results of the entire large sea surface were derived through extraction and synthesis. Figure 8 compares the backscattering coefficient results between the statistical SSA (theoretical results of SSA) and the proposed approach for both VV and HH polarizations with the incident azimuth of 0 • (upwind direction). The results of the proposed approach were averaged over 30 samples. As shown in Figure 8, good agreement exists between the corresponding results for all the incident angles, which suggests a good backscattering coefficient estimation of the sea surface at Ku-band with a wind speed of 5 m/s when using the proposed approach.

Numerical Simulations and Discussion
In this section, some numerical simulations are illustrated for sea surface scattering by applying the proposed method, including back and bi-static scattering coefficients varying with incident or scattering angles. Moreover, in order to demonstrate the validity of the proposed method it is compared to statistical SSA.
First, the backscattering simulation of sea surface was carried out at the Ku-band (14.0 GHz) with a wind speed of 5 m/s, which was chosen for the comparison with the statistical SSA results reported previously [27]. In this simulation, the large mesh grids were set to 0.1 × 0.1 m 2 , and the small grids were 0.001 × 0.001 m 2 to accurately estimate the Bragg scattering contribution. During the scattering field dataset construction, the typical cell included 100 large grid cells of the all possible inclinations for EM scattering calculation. In the implementation process, there were 1024 sampling points in both x and y directions, and the simulation area of the sea surface was 102.4 × 102.4 m 2 . The Elfouhaily model [28] was adopted for the sea roughness spectrum, which was developed based on available field and wave-tank measurements. The Elfouhaily model is supported by strong physical arguments, contrary to other spectra that are mostly empirical. The relative dielectric constant of the sea water was calculated using the Klein dielectric constant model [29] at 20 °C and 35 ‰ salinity.
With a 5 m/s wind speed, the variation ranges of the normal vector inclination angles, θ and φ versus the x-and y-axis, were 61-117 and 67-113°, respectively. After the dataset of the typical cell of the all possible inclinations was established for EM scattering calculation, interpolation was performed to obtain a much more delicate variation description with the interval of 0.1°. Finally, the scattering results of the entire large sea surface were derived through extraction and synthesis. Figure 8 compares the backscattering coefficient results between the statistical SSA (theoretical results of SSA) and the proposed approach for both VV and HH polarizations with the incident azimuth of 0° (upwind direction). The results of the proposed approach were averaged over 30 samples. As shown in Figure 8, good agreement exists between the corresponding results for all the incident angles, which suggests a good backscattering coefficient estimation of the sea surface at Ku-band with a wind speed of 5 m/s when using the proposed approach.
For a further verification, the backscattering results obtained using the two methods with a 15 m/s wind speed were compared (Figure 9). In this simulation, the other parameters and the operations remained the same with the above comparison, except for the wind speed-related values, such as the variation in the ranges of the normal vector inclination angles. With a 15 m/s wind speed, range in values changed to 56-124 a 61-121° for θ and φ , respectively. Again, the results of the proposed approach were averaged over 30 samples. From the comparison, the results obtained from both methods matched well for all incident angles and both VV and HH polarizations. This good agreement indicates the good performance of the proposed method for the EM scattering estimation of the sea surface under different sea conditions.  The above results only involve the backscattering calculations. The following simulations expand the sea surface scattering estimation to bistatic scattering to further validate the proposed method. First, a general forward-backward configuration of the bistatic scattering was considered, where the z-axis, the incident wave vectors, and the scattered wave vectors were in the same plane. The bistatic scattering simulation of sea surface was carried out at the Ku-band (14.0 GHz) with a wind speed of 5 m/s. Parameters such as sea spectrum, dielectric constant, and sampling intervals, were kept the same. The other parameters were fixed as follows: the transmitter incident and azimuth angles were 50° and 0°, respectively; the azimuth relative to wind direction was equal to 0 (upwind case); the receiver azimuth was set to180°; and the receiver incident angle s θ varied from −90° to 90°. Figure 10 compares the forward-backward bistatic scattering results between the proposed method and the statistical SSA, which is derived from Awada et al. [12]. As observed in the figure, the maximum energy was received around the specular direction 50°, which is a logical result because this is the true specular direction as given by Snell's law. This maximum decreased when the wind speed increased. On the other hand, good agreement was observed between the bistatic scattering results from both methods for all incident angles and both HH and VV polarizations. This agreement suggests that the proposed approach has good estimation ability for sea surface bistatic scattering at Ku-band at a wind speed of 5 m/s. To further illustrate the performance of the proposed approach for bistatic scattering simulation, another bistatic scattering simulation configuration was examined. In this simulation, parameters were the same as for the above example, except that the wind speed was changed to 15 The above results only involve the backscattering calculations. The following simulations expand the sea surface scattering estimation to bistatic scattering to further validate the proposed method. First, a general forward-backward configuration of the bistatic scattering was considered, where the z-axis, the incident wave vectors, and the scattered wave vectors were in the same plane. The bistatic scattering simulation of sea surface was carried out at the Ku-band (14.0 GHz) with a wind speed of 5 m/s. Parameters such as sea spectrum, dielectric constant, and sampling intervals, were kept the same. The other parameters were fixed as follows: the transmitter incident and azimuth angles were 50 • and 0 • , respectively; the azimuth relative to wind direction was equal to 0 (upwind case); the receiver azimuth was set to180 • ; and the receiver incident angle θ s varied from −90 • to 90 • . Figure 10 compares the forward-backward bistatic scattering results between the proposed method and the statistical SSA, which is derived from Awada et al. [12]. As observed in the figure, the maximum energy was received around the specular direction 50 • , which is a logical result because this is the true specular direction as given by Snell's law. This maximum decreased when the wind speed increased. On the other hand, good agreement was observed between the bistatic scattering results from both methods for all incident angles and both HH and VV polarizations. This agreement suggests that the proposed approach has good estimation ability for sea surface bistatic scattering at Ku-band at a wind speed of 5 m/s.  The above results only involve the backscattering calculations. The following simulations expand the sea surface scattering estimation to bistatic scattering to further validate the proposed method. First, a general forward-backward configuration of the bistatic scattering was considered, where the z-axis, the incident wave vectors, and the scattered wave vectors were in the same plane. The bistatic scattering simulation of sea surface was carried out at the Ku-band (14.0 GHz) with a wind speed of 5 m/s. Parameters such as sea spectrum, dielectric constant, and sampling intervals, were kept the same. The other parameters were fixed as follows: the transmitter incident and azimuth angles were 50° and 0°, respectively; the azimuth relative to wind direction was equal to 0 (upwind case); the receiver azimuth was set to180°; and the receiver incident angle s θ varied from −90° to 90°. Figure 10 compares the forward-backward bistatic scattering results between the proposed method and the statistical SSA, which is derived from Awada et al. [12]. As observed in the figure, the maximum energy was received around the specular direction 50°, which is a logical result because this is the true specular direction as given by Snell's law. This maximum decreased when the wind speed increased. On the other hand, good agreement was observed between the bistatic scattering results from both methods for all incident angles and both HH and VV polarizations. This agreement suggests that the proposed approach has good estimation ability for sea surface bistatic scattering at Ku-band at a wind speed of 5 m/s. To further illustrate the performance of the proposed approach for bistatic scattering simulation, another bistatic scattering simulation configuration was examined. In this simulation, parameters were the same as for the above example, except that the wind speed was changed to 15 To further illustrate the performance of the proposed approach for bistatic scattering simulation, another bistatic scattering simulation configuration was examined. In this simulation, parameters were the same as for the above example, except that the wind speed was changed to 15 m/s, the transmitter incident and azimuth angles were 40 • and 0 • , respectively, the receiver incident was set to 40 • , and the receiver azimuth angle φ s varied from 0 • to 180 • . Figure 11 compares the average NRCS between statistical SSA and the proposed approach. From the comparison, the results obtained with both methods were similar for all azimuth angles and both VV and HH polarizations. This good performance again shows the effectiveness of the proposed approach in bistatic sea surface scattering calculation for different scattering angles and sea conditions. Remote Sens. 2018, 10, x FOR PEER REVIEW 11 of 14 m/s, the transmitter incident and azimuth angles were 40° and 0°, respectively, the receiver incident was set to 40°, and the receiver azimuth angle s φ varied from 0° to 180°. Figure 11 compares the average NRCS between statistical SSA and the proposed approach. From the comparison, the results obtained with both methods were similar for all azimuth angles and both VV and HH polarizations. This good performance again shows the effectiveness of the proposed approach in bistatic sea surface scattering calculation for different scattering angles and sea conditions. As wind direction is also an important factor in sea surface scattering, investigating the performance of the proposed method under different wind directions was necessary. Given the relative geometry relationship between the incident azimuth and wind direction, the variation in backscattering coefficients as a function of azimuth angle of incident wave was calculated. Figure 12 compares the backscattering results as a function of azimuth angle of incident wave between the statistical SSA and the proposed method, where the incident wave was 14.0 GHz, the wind speed was 5 m/s, and the incident pitch angle was 60°. From the comparison, we observed that backscattering coefficients change with wind directions in a range of about 4 dB, and the results obtained with both methods were well-matched. Comparisons of the cross-polarization bistatic scattering coefficients between the statistical SSA (data derived from Awada et al. [12]) and the proposed method were also completed ( Figure  13). These comparisons were carried out in two different cases: (1) the transmitter incident and azimuth angles were 40° and 0°, respectively; the receiver azimuth was set to 45°; and the receiver incident angle s θ varied from 0° to 90°; and (2) the transmitter incident and azimuth angles were 40° and 0°, respectively; the receiver incident angle was set to 40°; and the receiver azimuth s φ varied from 0° to180°. The cross-polarization bistatic scattering coefficients derived from the two Figure 11. Comparison of the bistatic scattering results versus azimuth angles from the statistical SSA and the proposed method at a wind speed of 15 m/s. As wind direction is also an important factor in sea surface scattering, investigating the performance of the proposed method under different wind directions was necessary. Given the relative geometry relationship between the incident azimuth and wind direction, the variation in backscattering coefficients as a function of azimuth angle of incident wave was calculated. Figure 12 compares the backscattering results as a function of azimuth angle of incident wave between the statistical SSA and the proposed method, where the incident wave was 14.0 GHz, the wind speed was 5 m/s, and the incident pitch angle was 60 • . From the comparison, we observed that backscattering coefficients change with wind directions in a range of about 4 dB, and the results obtained with both methods were well-matched.
Remote Sens. 2018, 10, x FOR PEER REVIEW 11 of 14 m/s, the transmitter incident and azimuth angles were 40° and 0°, respectively, the receiver incident was set to 40°, and the receiver azimuth angle s φ varied from 0° to 180°. Figure 11 compares the average NRCS between statistical SSA and the proposed approach. From the comparison, the results obtained with both methods were similar for all azimuth angles and both VV and HH polarizations. This good performance again shows the effectiveness of the proposed approach in bistatic sea surface scattering calculation for different scattering angles and sea conditions. As wind direction is also an important factor in sea surface scattering, investigating the performance of the proposed method under different wind directions was necessary. Given the relative geometry relationship between the incident azimuth and wind direction, the variation in backscattering coefficients as a function of azimuth angle of incident wave was calculated. Figure 12 compares the backscattering results as a function of azimuth angle of incident wave between the statistical SSA and the proposed method, where the incident wave was 14.0 GHz, the wind speed was 5 m/s, and the incident pitch angle was 60°. From the comparison, we observed that backscattering coefficients change with wind directions in a range of about 4 dB, and the results obtained with both methods were well-matched. Comparisons of the cross-polarization bistatic scattering coefficients between the statistical SSA (data derived from Awada et al. [12]) and the proposed method were also completed ( Figure  13). These comparisons were carried out in two different cases: (1) the transmitter incident and azimuth angles were 40° and 0°, respectively; the receiver azimuth was set to 45°; and the receiver incident angle s θ varied from 0° to 90°; and (2) the transmitter incident and azimuth angles were 40° and 0°, respectively; the receiver incident angle was set to 40°; and the receiver azimuth s φ varied from 0° to180°. The cross-polarization bistatic scattering coefficients derived from the two Comparisons of the cross-polarization bistatic scattering coefficients between the statistical SSA (data derived from Awada et al. [12]) and the proposed method were also completed ( Figure 13).
These comparisons were carried out in two different cases: (1) the transmitter incident and azimuth angles were 40 • and 0 • , respectively; the receiver azimuth was set to 45 • ; and the receiver incident angle θ s varied from 0 • to 90 • ; and (2) the transmitter incident and azimuth angles were 40 • and 0 • , respectively; the receiver incident angle was set to 40 • ; and the receiver azimuth φ s varied from 0 • to180 • . The cross-polarization bistatic scattering coefficients derived from the two approaches were also well-aligned. Moreover, in Figure 13b, as θ i and θ s had the same value, the curves for VH and HV polarizations coincide due to the theoretical reason of SSA.
Remote Sens. 2018, 10, x FOR PEER REVIEW 12 of 14 approaches were also well-aligned. Moreover, in Figure 13b, as i θ and s θ had the same value, the curves for VH and HV polarizations coincide due to the theoretical reason of SSA. Through the above scattering result comparisons, we demonstrated that the proposed novel realization of SSA has the same accuracy for sea surface scattering for backscattering, bistatic scattering, and both polarizations. In addition to the accuracy attained, the computation efficiency of the scattering calculation for a large sea surface was significantly improved due to the combination of the dataset and field synthesis. In the above simulations, a total of 29 × 24 and 34 × 30 inclination situations were used for the 5 m/s and 15 m/s wind speeds, respectively, with an angle sampling interval of 2°. For each situation a 1 × 1 m 2 square area was calculated by the SSA. Accordingly, the total computation amount was equivalent to the EM scattering computation amount for the sea surface with an area of 29 × 24 m 2 and 34 × 30 m 2 , respectively. For a 100 × 100 m 2 area sea surface scattering simulation, the computation amount decreased to 7-10% of the traditional SSA realization. Importantly, this total computation amount was almost independent of the scale of the maritime scene to be simulated, and the performance of the efficiency improvement increased with the scale of the simulated sea surface. On the other hand, as the proposed approach uses the cooperation between large and small mesh grids, the one-eighth of wavelength sampling is no longer required and the computation memory is then greatly reduced, especially for high microwave bands and large sea surfaces. These improvements give the proposed method an outstanding advantage in large sea surface EM scattering simulation applications. Notably, the proposed method in this paper is only an efficiency-and memory-related improvement of the traditional numerical simulation of SSA. That is to say, this novel realization has the same accuracy and methodology limits as first-order SSA. The related simulations using SSA have drawbacks when tackling sea surface scattering under grazing angles and high sea states. In these relatively extreme cases, some extra effects should be considered, such as spiky and breaking waves, multiple scattering and the sheltering effect. These are all challenging problems to be tackled in the following work.

Conclusions
In this paper, a novel approach for SSA is proposed to significantly decrease the computation and computer memory requirements for sea surface scattering simulation. For this realization, the sea surface was decomposed into two scales, and each scale had its own spatial sample interval. Then, the scattering of the meshed facets of the cell surface was calculated under all the possible inclinations and a corresponding database was established. Once the dataset of the typical cell for EM scattering calculation of the all possible inclinations was established, the scattering results of the entire large sea surface were derived through extraction and synthesis operations. In the numerical Through the above scattering result comparisons, we demonstrated that the proposed novel realization of SSA has the same accuracy for sea surface scattering for backscattering, bistatic scattering, and both polarizations. In addition to the accuracy attained, the computation efficiency of the scattering calculation for a large sea surface was significantly improved due to the combination of the dataset and field synthesis. In the above simulations, a total of 29 × 24 and 34 × 30 inclination situations were used for the 5 m/s and 15 m/s wind speeds, respectively, with an angle sampling interval of 2 • . For each situation a 1 × 1 m 2 square area was calculated by the SSA. Accordingly, the total computation amount was equivalent to the EM scattering computation amount for the sea surface with an area of 29 × 24 m 2 and 34 × 30 m 2 , respectively. For a 100 × 100 m 2 area sea surface scattering simulation, the computation amount decreased to 7-10% of the traditional SSA realization. Importantly, this total computation amount was almost independent of the scale of the maritime scene to be simulated, and the performance of the efficiency improvement increased with the scale of the simulated sea surface. On the other hand, as the proposed approach uses the cooperation between large and small mesh grids, the one-eighth of wavelength sampling is no longer required and the computation memory is then greatly reduced, especially for high microwave bands and large sea surfaces. These improvements give the proposed method an outstanding advantage in large sea surface EM scattering simulation applications. Notably, the proposed method in this paper is only an efficiency-and memory-related improvement of the traditional numerical simulation of SSA. That is to say, this novel realization has the same accuracy and methodology limits as first-order SSA. The related simulations using SSA have drawbacks when tackling sea surface scattering under grazing angles and high sea states. In these relatively extreme cases, some extra effects should be considered, such as spiky and breaking waves, multiple scattering and the sheltering effect. These are all challenging problems to be tackled in the following work.

Conclusions
In this paper, a novel approach for SSA is proposed to significantly decrease the computation and computer memory requirements for sea surface scattering simulation. For this realization, the sea surface was decomposed into two scales, and each scale had its own spatial sample interval.
Then, the scattering of the meshed facets of the cell surface was calculated under all the possible inclinations and a corresponding database was established. Once the dataset of the typical cell for EM scattering calculation of the all possible inclinations was established, the scattering results of the entire large sea surface were derived through extraction and synthesis operations. In the numerical simulations, backscattering and bistatic scattering of the sea surface were simulated for different wind speeds, incident and scattering angles, and polarizations. From the comparisons with statistical SSA results, we demonstrated that the proposed approach possesses the same accuracy for relative sea surface scattering applications. In addition to the accuracy, the computation efficiency of the scattering calculation for a large sea surface was significantly improved, due to the combination of a database and field synthesis. Moreover, one-eighth of the wavelength sampling is no longer needed and the computation memory is then greatly reduced. With these improvements, this proposed approach has prominent advantages in sea surface EM scattering simulation applications, especially for high microwave bands and large sea surfaces. Unlike statistical SSA, which is an analytical method and can only be applied to the prediction of mean scattering coefficients from sea surfaces, the proposed method is an improved Monte-Carlo simulation of sea surface scattering, which includes sea surface generation and the EM scattering calculation processes from every surface section. This model can not only be used for mean scattering coefficients prediction, but can also be applied to relative simulations about the specific sea surface profile, such as composite scattering problems of sea surfaces with targets, synthetic aperture radar (SAR) image simulation of maritime scenes, and the dynamic scattering characteristics of the sea surface.