Efficient Numerical Full-Polarized Facet-Based Model for EM Scattering from Rough Sea Surface within a Wide Frequency Range

Abstract: A full-polarized facet based scattering model (FPFSM) for investigating the electromagnetic (EM) scattering by two-dimensional electrically large sea surfaces with high efficiency at high microwave bands is proposed. For this method, the scattering field over a large sea facet in a diffuse scattering region is numerically deduced according to the Bragg scattering mechanism. In regard to near specular directions, a novel approach is proposed to calculate the scattered field from a sea surface based on the second order small slope approximation (SSA-II), which saves computer memory considerably and is able to analyze the EM scattering by electrically large sea surfaces. The feasibility of this method in evaluating the radar returns from the sea surface is proved by comparing the normalized radar cross sections (NRCS) and the Doppler spectrum with the SSA-II. Then NRCS results in monostatic and bistatic configurations under different polarization states, scattering angles and wind speeds are analyzed as well as the Doppler spectrum at Ka-band. Numerical results show that the FPFSM is a reliable and efficient method to analyze the full-polarized scattering characteristics from electrically large sea surface within a wide frequency range.


Introduction
Full-polarized echoes from a sea's surface include more abundant information than single-polarized ones [1][2][3][4].Thus, the study of the full-polarized scattering characteristics from sea surfaces is of great value and will be of benefit in the detection and recognition of targets in a maritime environment.
The electromagnetic (EM) scattering by a rough sea surface has been carefully studied for a long time [5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22], and many methods have been proposed.Generally, due to the electrically large size of the sea surface at microwave bands, analytical approximation methods are widely employed, such as the Kirchhoff approximation (KA), the small perturbation method (SPM) [8], the two-scale model (TSM) [6] and the small-slope approximation (SSA) [9].All of these methods have some limitations.The KA is only valid for sea surfaces with small values for the ratio of wavelength to curvature radius.Based on KA, a faceted approach is proposed in Reference [17] to simulate the bistatic returns from the ocean surface which works well in a quasi-specular regime.As for the SPM, it is only valid for surfaces with small values for the ratio of roughness to the incident wavelength.Based on the theory of composite surface, the famous TSM is presented to deal with the scattering by rough sea surfaces.However, the TSM mainly suffers from the problem of how to distinguish between gravity waves and capillary waves.Based on the idea of TSM, several extended methods [18][19][20][21] have been proposed according to the Bragg scattering mechanism.In References [18,19], the Bragg-scattering is used to explain most features of sea surface backscatter at moderate incident angles.The non-Bragg component of backscatter is also considered at higher incident angles.In Franceschetti et al. [20], the synthetic aperture radar (SAR) echoes from the sea surface at moderate incidence angles are simulated based on the Bragg-scattering mechanism.Then a facet-based method is derived in Li et al. [21], which works well at small to moderate incident angles because different scattering mechanisms are considered.However, if only Bragg scattering is used, it can only be applied at moderate incident angles.Unlike the TSM, the SSA is a unifying theory that is able to deal with small roughness as well as large curvatures.Moreover, it has a wider application range [9].However, the SSA method requires that the sampling interval of rough surfaces must be less than one-eighth of the incident wavelength.Thus, it suffers from the problem of heavy consumption of computer memory.Based on the SSA, some efforts have been made to reduce the computational complexity.In Guérin et al. [23], the high-frequency limit of the second order SSA (SSA-II) is proposed, and it is able to calculate the co-polarized scattering coefficient at high microwave bands efficiently in both monostatic and bistatic configurations.Then in Guérin et al. [1], a simplified formulation for rough surface cross-polarized backscattering radar cross sections (RCS) was proposed, which makes a significant contribution.In Voronovich et al. [24], a theoretical model to predicate the scattering from sea surface at Ku and C bands is proposed based on the sea spectrum, which has good accuracy when compared with measured data.Then the simplified forms for the first three terms in the SSA series for incoherent scattering from a Gaussian correlation function surface were presented in Gilbert et al. [25].All these studies mainly focus on the prediction of the scattering coefficient from rough surfaces.
Nowadays, with the development of high-resolution radar, the averaged RCS results are not enough to address remote sensing problems.Faced with this problem, this paper proposed the full-polarized facet based scattering model (FPFSM), which is not only able to give the averaged RCS results but also able to give the full-polarized plural scattered field including amplitude and phase information.As for this method, the scattered field from a large sea facet in diffuse scattering directions was numerically calculated according to the Bragg scattering mechanism.Then in the near specular direction, a different way based on SSA-II is proposed to calculate the scattered field, which saves memory consumption to a large extent and makes it possible to analyze the full-polarized scattering by electrically large sea surfaces.In conclusion, the main contribution is that we proposed a facet-based method to investigate the EM scattering characteristic of electrically large sea surfaces.This method works well from L-band to Ka-band and also has good accuracy, much higher computational efficiency, and much lower memory consumption compared with the original SSA-II method.In addition, it is able to give the plural scattering field information (i.e., containing both amplitude and phase information) of each facet, which means it can be further applied to analyze the scattering from targets over sea surfaces and to the high-resolution of SAR imaging simulators.

Near Specular Reflection Direction
In this section, the method to predict the scattered field in near specular reflection directions is proposed based on the SSA-II.Assuming that the incident wave vector and the scattering wave vector are respectively k i = k ix x + k iy ŷ − q 0 ẑ and k s = k sx x + k sy ŷ + q 1 ẑ, then the scattering amplitude (SA) is [9] S PQ (k, k 0 ) = 2 √ q 1 q 0 q 1 +q 0 dr 4π 2 exp[−j(k − k 0 ) • r + j(q 1 + q 0 )h(r)] × B PQ (k, k 0 ) − j 4 M PQ (k, k 0 ; ξ)h(ξ) exp(jξ • r)dξ (1) where PQ= HH, VV, HV, VH represents the polarization mode, k 0 = k ix x + k iy ŷ and k = k sx x + k sy ŷ are the horizontal projections of the incident wave vector k i and the scattered wave vector k s , respectively; ξ = (ξ x , ξ y ), ξ x , ξ y are the wavenumbers along x and y directions; h(ξ) is the Fourier transform of the sea surface height h(r) with the expression where a(ξ) is the complex amplitude of the wave, which is related to the ocean spectrum; ω ξ is the gravity-capillarity dispersion relationship.The expressions for B(k, k 0 ) at different polarization modes are [9,24] where ε is the dielectric constant of sea water, which can be calculated with the method in Klein et al. [26]; are the vertical components of the incident waves in the air and sea water, respectively and q 2 = εk 2 i − k 2 0 , q 02 = εk 2 i − k 2 are the vertical components of the scattered waves in the air and the sea water, respectively; k = |k|, k 0 = |k 0 |, N = (0, 0, 1) is the unit vertical vector.Then, the expression of the expressions for B 2,PQ are [9,24] where Note that the imaginary parts of q 1 , q 01 q 2 , q 02 , q ξ1 , q ξ2 are all non-negative.
In fact, when generating a sea surface with M × N intervals, h(ξ) should be generated first, then the two-dimensional (2D) inverse fast Fourier transformation (IFFT) is applied to generate the sea height.For the SSA-II, the whole sea surface of interest is generally generated with an interval less than one-eighth of the incident wavelength, which requires very large computer memory availability at higher microwave bands.On the other hand, the second order scattering term M(k, k 0 ; ξ)h(ξ) exp(jξ • r)dξ enhances the computational complexity even further.
As for the near specular direction, large-scale waves are the main contribution to the scattering field [27][28][29].Thus, an approximated h (ξ) is used to replace h(ξ), and the relationship between them is where k cut is the cut-off wave number.Based on the analyses above, when analysing the scattering from a sea surface in near nadir direction, the surface can be generated by applying 2D IFFT to h (ξ).
In this process, according to Equation ( 12), the 2D matrix h (ξ) has the following form, i.e., where A is the matrix that contains the non-zero elements of h (ξ) with the size of M 1 × N 1 .Then the 2D IFFT of h (ξ) can be obtained by Equation ( 14) where B = F row 0 0 A 0 0 ; F 2 , F row , F col represent the 2D IFFT, the 1D IFFT along row and column directions, respectively.The consumed memory for the F row operation is O(M 1 • N).Then for the F col operation, the consumed memory can be as small as O(M).Thus, the proposed method decreases the consumed memory significantly compared with the original 2D IFFT operation, of which the consumed memory is O(M • N).That is to say, the proposed method has the ability to generate the electrically large sea surface.Then in regards to the calculation of the second order scattered field in Equation ( 2), a similar method can be used to replace the 2D IFFT operation.
With the method above, the normalized radar cross sections (NRCS) results in near specular directions are calculated by the proposed method with different cut-off wavenumbers k cut = k/3, k/4, k/5 and are compared with that calculated by the original SSA-II in Figure 1.Here, k is the incident wavenumber.In this simulation, M = N = 1024, the wind speed is 5 m/s, the incident frequency is 1 GHz.One can observe that the error increases with incident angles, i.e., the contribution to the total scattered field from small-scale waves increases with a larger incident angle.In addition, good agreement is shown within a wider numerical range of incident angles when a larger cut-off wavenumber is employed.If we imagine that it is infinite, the simplification effect will disappear and the method will be the same as the original SSA-II.Moreover, by comparing Figure 1a-c, it is observed that with the same cut-off wavenumber, good agreement is shown within a wider numerical range of incident angles for co-polarized modes.Thus, we can conclude that for a fixed incident angle, the cross-polarized scattered field comes from many more wave components than the co-polarized scattered field.Moreover, from these comparisons, we can observe that k cut = k/3 is appropriate because it works well in both co-polarization and cross-polarization modes in near specular directions.
In fact, a larger cut-off wavenumber will work better; however, the consumed memory increases simultaneously.Thus, considering the calculation accuracy and memory consumption, k cut = k/3 is selected.
disappear and the method will be the same as the original SSA-II.Moreover, by comparing Figure 1a-c, it is observed that with the same cut-off wavenumber, good agreement is shown within a wider numerical range of incident angles for co-polarized modes.Thus, we can conclude that for a fixed incident angle, the cross-polarized scattered field comes from many more wave components than the co-polarized scattered field.Moreover, from these comparisons, we can observe that cut / 3 k k  is appropriate because it works well in both co-polarization and cross-polarization modes in near specular directions.In fact, a larger cut-off wavenumber will work better; however, the consumed memory increases simultaneously.Thus, considering the calculation accuracy and memory consumption, cut / 3 k k  is selected.Then in Figure 2, the bistatic NRCS results calculated by the proposed method in near specular directions are compared with that calculated by the original SSA-II as well.The sampling numbers, the wind speed and direction, and the radar parameters are the same as Figure 1.The cut-off wavenumber is selected as / 3 k .The incident angle is i 40 From the comparison result in Figures 1 and 2, it is obvious that the method shown in this section has a good agreement with the original SSA-II in near specular directions.This good agreement also proves that even for cross-polarizations, the scattered field from the large-scale waves is the main contribution to the total field in near specular directions.Moreover, considering the method in this subsection, it significantly saves computer memory compared with the original SSA-II.Thus, the Then in Figure 2, the bistatic NRCS results calculated by the proposed method in near specular directions are compared with that calculated by the original SSA-II as well.The sampling numbers, the wind speed and direction, and the radar parameters are the same as Figure 1.The cut-off wavenumber is selected as k/3.The incident angle is θ i = 40 • , the azimuth angle is ϕ i = 0 • , the scattering azimuth angle is ϕ s = 0 • and the scattering angle θ s varies from −80 • to 80 • , and particularly, θ s = θ i is the specular direction and θ s = −θ i is the backward scattering direction.The definition of the angles is shown in Figure 3, and for such a coordinate, the relation between ki , ks and the four angles are ki = x sin θ i cos φ i + ŷ sin θ i sin φ i − ẑ cos θ i ks = x sin θ s cos φ s + ŷ sin θ s sin φ s + ẑ cos θ s (15) Remote Sens. 2019, 11, x FOR PEER REVIEW 6 of 18 proposed method makes it possible to analyze the scattered field from very electrically-large rough sea surfaces on personal computers.From the comparison result in Figures 1 and 2, it is obvious that the method shown in this section has a good agreement with the original SSA-II in near specular directions.This good agreement also proves that even for cross-polarizations, the scattered field from the large-scale waves is the main contribution to the total field in near specular directions.Moreover, considering the method in this subsection, it significantly saves computer memory compared with the original SSA-II.Thus, the proposed method makes it possible to analyze the scattered field from very electrically-large rough sea surfaces on personal computers.
proposed method makes it possible to analyze the scattered field from very electrically-large rough sea surfaces on personal computers.

Diffuse Scattering Region
In the off-nadir region, first order Bragg scattering is the main contribution to the co-polarized scattered field [18][19][20].However, both the tilt modulation and the second order Bragg scattering contribute to the depolarization.Even though the second order Bragg scattering mechanism contributed to the cross-polarization scattering field, there was only about a 2~4 dB difference if it is neglected when 10 15m/s u  [2].Thus, for the sake of computational efficiency, the two-scale model is modified with the first order Bragg scattering mechanism and the tilt modulation effect.In this way, both the effects of tilting modulation and local small-scale capillary waves are taken into consideration, which is reasonable [30].
In the following, the right-handed global and local coordinate systems are adopted and shown in Figure 3.For global one, the z-axis points vertically upward and the ˆxOy plane lies on the mean sea level plane.The local coordinate system is ˆˆ-   O x y n.  O and n are the center and the normal vector of a sea facet, respectively.Then based on Fuks' formula, the scattered field from an arbitrary facet can be written as [31,32]   where i E is the incident wave, ( )  r is the small-scale capillary wave riding on the tilt facet at position r ; PQ   is the polarization factors which is related to the relative dielectric constant  of sea water and Fresnel reflection coefficients; are the polarization of incident and scattered waves, respectively;  q is the projection of the k k on the tilted plane, and k is the wave number of electromagnetic wave.
It was generally accepted that the main contribution of these types of short waves (tilted by the large-scale profiles) to the microwave radar return depends on the Bragg resonant mechanism between the EM waves and sea surface components of Bragg wavelength [33,34].Thus, it is

Diffuse Scattering Region
In the off-nadir region, first order Bragg scattering is the main contribution to the co-polarized scattered field [18][19][20].However, both the tilt modulation and the second order Bragg scattering contribute to the depolarization.Even though the second order Bragg scattering mechanism contributed to the cross-polarization scattering field, there was only about a 2~4 dB difference if it is neglected when u 10 < 15 m/s [2].Thus, for the sake of computational efficiency, the two-scale model is modified with the first order Bragg scattering mechanism and the tilt modulation effect.In this way, both the effects of tilting modulation and local small-scale capillary waves are taken into consideration, which is reasonable [30].
In the following, the right-handed global and local coordinate systems are adopted and shown in Figure 3.For global one, the z-axis points vertically upward and the xO ŷ plane lies on the mean sea level plane.The local coordinate system is O − x ŷ n.O and n are the center and the normal vector of a sea facet, respectively.Then based on Fuks' formula, the scattered field from an arbitrary facet can be written as [31,32] where E i is the incident wave, ζ(r) is the small-scale capillary wave riding on the tilt facet at position r; F PQ is the polarization factors which is related to the relative dielectric constant ε of sea water and Fresnel reflection coefficients; P = H, V and Q = H, V are the polarization of incident and scattered waves, respectively; q ⊥ is the projection of the q = k(k s − k i ) on the tilted plane, and k is the wave number of electromagnetic wave.It was generally accepted that the main contribution of these types of short waves (tilted by the large-scale profiles) to the microwave radar return depends on the Bragg resonant mechanism between the EM waves and sea surface components of Bragg wavelength [33,34].Thus, it is reasonable to suppose that only the ripples traveling along the line of the incident-scatter vector q are visible to the radar sensor.Then, the high-frequency microstructure of each facet can be represented by a single cosinoidal ripple that has the Bragg resonant wavelength, as where r 0 is the center of the large facet, B(κ) = 2πS(κ) represents the amplitude of the short waves, wherein S(κ) is the anisotropy capillary spectrum located in the higher part of the sea spectrum.
The direction of κ is the same with q ⊥ and its amplitude can be obtained based on the Bragg wavelength λ b = λ/sin θ q , wherein λ is the incident wavelength, θ q is the angle between q and the normal vector n of the facet.In addition, in Figure 3, κ x and κ y are the projections of the Bragg wave vector κ along the x and ŷ directions, respectively.Then the integration term in Equation ( 16) over a tilt facet is where (x , y ) is the local coordinate of a point r on the tilt facet referring to its center; q 1 , q 2 are the projections of q ⊥ along x and ŷ directions, respectively; δx = 1 + z 2 x ∆x, δy = 1 + z 2 y ∆y, ∆x and ∆y are the sampling intervals of the sea surface along x and y directions in the xo ŷ plane, respectively, and z x and z y are the slopes of the tilt facet along x and y directions, respectively.The integration in Equation ( 18) can be further simplified as Up to now, the scattered field from a facet is numerically derived based on the Bragg scattering mechanism and it can be calculated with Equations ( 16), (18), and (19), which depend on the slopes of the facet and the wind speed, wind direction, and the radar configuration.In addition, the capillary wave tilted by the large-scale gravity waves are assumed as cosinoidal ripples with specific Bragg resonant wavelength, thus, the sea surface can be generated with an interval of ∆x = ∆y = γλ(γ > 1), which greatly decreases memory consumption and enhances the computational efficiency.It needs to be noted that this method can only be employed in diffuse scattering directions.

The Complete Full-Polarized Facet Based Scattering Model
Then for a sea surface generated with an interval of λ/8, which is mentioned as sea-1, the total scattered field can be calculated with the method in Section 2.1 by summing the vector field from different small facets.Then, we can also obtain a sea surface with an interval of γλ(γ > 1) on the basis of sea-1, and this one is named as sea-2.Assuming two matrixes η 1 and η 2 represent the height of sea-1 and sea-2, respectively.Then the relationship between them is where γ is an integer, M, N are the sampling points of sea-1 along x and y directions.Then the scattered field from sea-2 can be calculated with the method in Section 2.2.Up to now, the total field from a sea surface can be obtained by the summation of the fields from all the facets on the sea surface, Remote Sens. 2019, 11, 75 Finally, the NRCS σ FPFSM can be obtained by where A s is the area of the sea surface.In fact, E sea-1 can be regarded as 0 in diffuse specular directions, which means that in most radar configurations, E sea-1 is no longer needs to be calculated and computational efficiency will be greatly enhanced by only calculating E sea-2 .

Performance of Full-Polarized Facet Based Scattering Model
In this section, the validity of the proposed FPFSM is verified by comparing the NRCS results calculated by the proposed method and the original SSA-II for different radar configurations, incident frequencies, and polarizations.It should be noted that for low grazing angles, the shading effect and the existence of irregular waves (such as breaking waves) have an effect on the NRCS results.The proposed method does not consider these effects because the main purpose of this paper is not to solve the problems at low grazing angles but to propose a facet-based method to investigate the full-polarized EM scattering characteristic of electrically large sea surfaces efficiently.Therefore, the numerical range of the incident angle is 0 • to 80 • in monostatic configuration, and the scattering angle varies from −80 • to 80 • in bistatic configuration.In addition, the calculation time and the memory consumption are also compared to show the higher efficiency and much lower memory consumption of FPFSM.Furthermore, the Doppler spectrum results are also shown to further demonstrate its good performance.

Normalized Radar Cross Section Results for Different Wave Bands
In this paper, the two-dimensional sea surfaces are all generated using the Elfouhaily spectrum [35] accompanied by Longuet-Higgins' directional function [36].First, Figure 4a-c shows the averaged NRCS results calculated by FPFSM with different samples.The tapered incident wave [37] is employed, and the tapering parameter is selected as L/4 with L the dimension of the sea surface, the incident frequency is 1 GHz.The size of the sea surface is 40 m × 40 m, the wind speed is u 10 = 5 m/s at 10 m height above the mean sea level in the upwind direction, and the Elfouhaily spectrum as a function of wavenumber for such a wind speed is shown in Figure 4d.It can be observed that the NRCS results from one sample show obvious differences with the averaged results, and the averaged NRCS results tend to be stable for a fixed incident angle with the increasing samples.Thus, the averaged NRCS results from 50 samples are used in this paper to do the comparison with the SSA-II method.
Finally, the NRCS FPFSM  can be obtained by (24) where s A is the area of the sea surface.In fact, sea-1 E can be regarded as 0 in diffuse specular directions, which means that in most radar configurations, sea-1 E is no longer needs to be calculated and computational efficiency will be greatly enhanced by only calculating sea-2 E .

Performance of Full-polarized Facet based Scattering Model
In this section, the validity of the proposed FPFSM is verified by comparing the NRCS results calculated by the proposed method and the original SSA-II for different radar configurations, incident frequencies, and polarizations.It should be noted that for low grazing angles, the shading effect and the existence of irregular waves (such as breaking waves) have an effect on the NRCS results.The proposed method does not consider these effects because the main purpose of this paper is not to solve the problems at low grazing angles but to propose a facet-based method to investigate the full-polarized EM scattering characteristic of electrically large sea surfaces efficiently.Therefore, the numerical range of the incident angle is 0° to 80° in monostatic configuration, and the scattering angle varies from −80° to 80° in bistatic configuration.In addition, the calculation time and the memory consumption are also compared to show the higher efficiency and much lower memory consumption of FPFSM.Furthermore, the Doppler spectrum results are also shown to further demonstrate its good performance.

Normalized Radar Cross Section Results for Different Wave Bands
In this paper, the two-dimensional sea surfaces are all generated using the Elfouhaily spectrum [35] accompanied by Longuet-Higgins' directional function [36].First, Figure 4a-c shows the averaged NRCS results calculated by FPFSM with different samples.The tapered incident wave [37] is employed, and the tapering parameter is selected as / 4 L with L the dimension of the sea surface, the incident frequency is 1 GHz.The size of the sea surface is 40 m 40 m  , the wind speed is 10 5 m/s u  at 10 m height above the mean sea level in the upwind direction, and the Elfouhaily spectrum as a function of wavenumber for such a wind speed is shown in Figure 4d.It can be observed that the NRCS results from one sample show obvious differences with the averaged results, and the averaged NRCS results tend to be stable for a fixed incident angle with the increasing samples.Thus, the averaged NRCS results from 50 samples are used in this paper to do the comparison with the SSA-II method.Then, Figures 5-7 illustrate the comparisons of NRCS results that calculated by FPFSM and SSA-II at ultra high frequency (UHF) band (300 MHz), L-band (1 GHz), and C-band (8 GHz) for monostatic and bistatic configurations, and the NRCS results are averaged over 50 samples of sea surfaces.The size of the sea surface is 40 m 40 m  , the wind speed is 10 5 m/s u  in upwind direction.Then in Figure 8, the NRCS results for higher wind speeds ( 10 8 m/s,10 m/s u  ) are calculated and compared with the SSA-II as well.The good agreement of the results calculated by the two methods indicates that the good performance of the proposed method.Moreover, the NRCS results as a function of incident angle at C band are calculated with FPFSM and SSA-II and are compared with real SAR data in Figure 6a.The real SAR data are obtained with the CMOD-5 model [38] for horizontal-transmit-horizontal-receive (HH) polarization mode.Then, one can calculate the vertical-transmit-vertical -receive (VV) polarized results based on the polarization ratio [39].As for cross-polarized case, the real SAR data are calculated with the crosspolarization geophysical model function for C-band radar backscattering [40].From the simulated results, we can observe that the co-polarized NRCS results calculated by the FPFSM and the SSA-II agree well with the measured data.However, a discrepancy occurs of about 2~5 dB for the crosspolarization mode.The obvious difference appears in other researches as well the reason for which is the measurement of cross-polarized radar return from the ocean surface is always a challenge due to the weakness of this component [41].
On the other hand, the computation time for calculating the NRCS results in Figures 5-7 with one sample of sea surface by the two methods with an angle interval of 1 are shown in Table 1.These results are calculated by central processing unit (CPU) under the configuration: Intel(R) Core(TM) i5-760@3.3GHz.The comparison of the computation time indicates higher computing efficiency of FPFSM compared with SSA-II.In addition, the comparison of the peak memory consumption is given in Table 2. Noting that the memory consumption is the same in monostatic and bistatic configurations for a fixed frequency and surface size, we only show the comparison in monostatic configuration.In addition, for the FPFSM, the detailed implementation methods are not the same for different scattering directions, thus, the consumed memory of FPFSM in the near specular reflection direction and diffuse scattering direction are shown separately.The comparison result indicates that the proposed method saves a significant memory when analysing the scattering characteristic from sea surfaces, which means it is able to be applied to super electrically large sea surfaces.
Moreover, the mean and the maximum errors between the FPFSM and the SSA-II are shown in Table 3 to quantify the goodness of fit between the two methods.The mean error and the maximum error are defined as Then, Figures 5-7 illustrate the comparisons of NRCS results that calculated by FPFSM and SSA-II at ultra high frequency (UHF) band (300 MHz), L-band (1 GHz), and C-band (8 GHz) for monostatic and bistatic configurations, and the NRCS results are averaged over 50 samples of sea surfaces.The size of the sea surface is 40 m × 40 m, the wind speed is u 10 = 5 m/s in upwind direction.Then in Figure 8, the NRCS results for higher wind speeds (u 10 = 8 m/s, 10 m/s) are calculated and compared with the SSA-II as well.The good agreement of the results calculated by the two methods indicates that the good performance of the proposed method.
Moreover, the NRCS results as a function of incident angle at C band are calculated with FPFSM and SSA-II and are compared with real SAR data in Figure 6a.The real SAR data are obtained with the CMOD-5 model [38] for horizontal-transmit-horizontal-receive (HH) polarization mode.Then, one can calculate the vertical-transmit-vertical -receive (VV) polarized results based on the polarization ratio [39].As for cross-polarized case, the real SAR data are calculated with the cross-polarization geophysical model function for C-band radar backscattering [40].From the simulated results, we can observe that the co-polarized NRCS results calculated by the FPFSM and the SSA-II agree well with the measured data.However, a discrepancy occurs of about 2~5 dB for the cross-polarization mode.
The obvious difference appears in other researches as well the reason for which is the measurement of cross-polarized radar return from the ocean surface is always a challenge due to the weakness of this component [41].
On the other hand, the computation time for calculating the NRCS results in Figures 5-7 with one sample of sea surface by the two methods with an angle interval of 1 • are shown in Table 1.These results are calculated by central processing unit (CPU) under the configuration: Intel(R) Core(TM) i5-760@3.3GHz.The comparison of the computation time indicates higher computing efficiency of FPFSM compared with SSA-II.In addition, the comparison of the peak memory consumption is given in Table 2. Noting that the memory consumption is the same in monostatic and bistatic configurations for a fixed frequency and surface size, we only show the comparison in monostatic configuration.In addition, for the FPFSM, the detailed implementation methods are not the same for different scattering directions, thus, the consumed memory of FPFSM in the near specular reflection direction and diffuse scattering direction are shown separately.The comparison result indicates that the proposed method saves a significant memory when analysing the scattering characteristic from sea surfaces, which means it is able to be applied to super electrically large sea surfaces.
Moreover, the mean and the maximum errors between the FPFSM and the SSA-II are shown in Table 3 to quantify the goodness of fit between the two methods.The mean error and the maximum error are defined as where σ SSA−II is the NRCS results calculated with SSA-II.It is shown that the maximum errors are about 2~3 dB and mean errors are only about 1 dB, which indicates the good performance of FPFSM.

Comparison of Doppler Spectrum
Apart from the NRCS results, the Doppler spectrum of the sea surface is also of great value to obtain the sea condition information.The Doppler effect is due to the change in frequency of a wave in relation to an observer who is moving relative to the wave source.In regards to the ocean surface, it is caused by the change of the sea surface velocity projected along the line-of-sight of the radar.As two parameters of the Doppler spectrum, the Doppler shift reflects the line-of-sight velocity of the scatterers weighted by their contribution to the backscattered power [42], and the bandwidth of the Doppler spectrum is determined by the variance of the velocity distribution of the scattering facets at the sea surface [43].In this subsection, the proposed method is used to analyze the Doppler spectrum of the time-varying sea surface at L-band with a standard spectral estimation technique [44,45].It can be given as

Comparison of Doppler Spectrum
Apart from the NRCS results, the Doppler spectrum of the sea surface is also of great value to obtain the sea condition information.The Doppler effect is due to the change in frequency of a wave in relation to an observer who is moving relative to the wave source.In regards to the ocean surface, it is caused by the change of the sea surface velocity projected along the line-of-sight of the radar.As two parameters of the Doppler spectrum, the Doppler shift reflects the line-of-sight velocity of the scatterers weighted by their contribution to the backscattered power [42], and the bandwidth of the Doppler spectrum is determined by the variance of the velocity distribution of the scattering facets at the sea surface [43].In this subsection, the proposed method is used to analyze the Doppler spectrum of the time-varying sea surface at L-band with a standard spectral estimation technique [44,45].It can be given as

Comparison of Doppler Spectrum
Apart from the NRCS results, the Doppler spectrum of the sea surface is also of great value to obtain the sea condition information.The Doppler effect is due to the change in frequency of a wave in relation to an observer who is moving relative to the wave source.In regards to the ocean surface, it is caused by the change of the sea surface velocity projected along the line-of-sight of the radar.As two parameters of the Doppler spectrum, the Doppler shift reflects the line-of-sight velocity of the scatterers weighted by their contribution to the backscattered power [42], and the bandwidth of the Doppler spectrum is determined by the variance of the velocity distribution of the scattering facets at the sea surface [43].In this subsection, the proposed method is used to analyze the Doppler spectrum of the time-varying sea surface at L-band with a standard spectral estimation technique [44,45].It can be given as where E s PQ (t) is the scattered field calculated with Equation ( 21) at time t, the angular bracket • means the average from different samples of sea surfaces, T is the total evolution time for a time-varying sea surface.
Figure 9 illustrates the normalized Doppler spectrum averaged over 100 samples of sea surfaces for different polarizations and incident angles.In the simulation, the wind speed is 5 m/s in the upwind direction, the time step is set as 0.025 s and the total evolution time is 3.2 s.The size of the sea surface is 50 m × 50 m.The frequency of the incident wave is 1 GHz.As indicated in Figure 8, one can observe the good performance of the proposed method on investigating the Doppler spectrum for different radar configurations and polarizations.

Analyses of Electromagnetic Scattering Characteristics at Ka-Band
In this section, the FPFSM is further applied to Ka-band (30 GHz).The NRCS results and the Doppler spectrums from sea surfaces at Ka-band are calculated and analyzed by the proposed method.

Normalized Radar Cross Section under Monostatic Configuration
Firstly, the backscattered NRCS results calculated by the FPFSM are given in Figure 10.The size of the sea surface and the wind speed is the same as that in Section 3.1.At Ka-band, the SSA-II is difficult to implement due to the very large requirement for computer memory.Thus, for the comparison, the co-polarized results calculated with FPFSM are compared with those generated by the first order small slope approximation (SSA-I) with the spectral decomposition method [46], and the cross-polarized results are compared with the data in Guérin et al. [1].From Figure 10, it can be observed that the two models are in good agreement, which indicates the good performance of FPFSM in evaluating the radar returns from the sea surface at higher microwave bands.

Analyses of Electromagnetic Scattering Characteristics at Ka-Band
In this section, the FPFSM is further applied to Ka-band (30 GHz).The NRCS results and the Doppler spectrums from sea surfaces at Ka-band are calculated and analyzed by the proposed method.

Normalized Radar Cross Section under Monostatic Configuration
Firstly, the backscattered NRCS results calculated by the FPFSM are given in Figure 10.The size of the sea surface and the wind speed is the same as that in Section 3.1.At Ka-band, the SSA-II is difficult to implement due to the very large requirement for computer memory.Thus, for the comparison, the co-polarized results calculated with FPFSM are compared with those generated by the first order small slope approximation (SSA-I) with the spectral decomposition method [46], and the cross-polarized results are compared with the data in Guérin et al. [1].From Figure 10, it can be observed that the two models are in good agreement, which indicates the good performance of FPFSM in evaluating the radar returns from the sea surface at higher microwave bands.
difficult to implement due to the very large requirement for computer memory.Thus, for the comparison, the co-polarized results calculated with FPFSM are compared with those generated by the first order small slope approximation (SSA-I) with the spectral decomposition method [46], and the cross-polarized results are compared with the data in Guérin et al. [1].From Figure 10, it can be observed that the two models are in good agreement, which indicates the good performance of FPFSM in evaluating the radar returns from the sea surface at higher microwave bands.

Normalized Radar Cross Section under Bistatic Configuration
Furthermore, the NRCS results at Ka-band under forward-backward configuration is shown in Figure 11a.Here, the incident angle is i 40    , and the scattering angle s  varies from -80 to 80 .
Specifically, in this configuration, one can observe that the maximum energy is received around the

Normalized Radar Cross Section under Bistatic Configuration
Furthermore, the NRCS results at Ka-band under forward-backward configuration is shown in Figure 11a.Here, the incident angle is θ i = 40 • , and the scattering angle θ s varies from −80 • to 80 • .Specifically, in this configuration, one can observe that the maximum energy is received around the specular direction for co-polarized cases, which is a logical result.However, for cross-polarized results, the maximum value does not appear in the specular direction.In addition, unlike the backscattering configuration, the NRCS results for vertical-transmit-horizontal-receive (HV) polarization are not exactly the same as the horizontal-transmit-vertical-receive (VH) polarized results except for the specular direction.
Furthermore, the NRCS results in fully bistatic configurations at Ka-band are investigated as well.Figure 11b,c gives the NRCS with the incident direction θ i = 40 o , ϕ i = 0 o and the scattering direction: (b) θ s = −80 0 ∼ 80 0 , ϕ s = 45 o ; (c) θ s = −80 0 ∼ 80 0 , ϕ s = 135 o .It is interesting to find that the discrepancy between the VH and HV polarizations became wider in these configurations.In addition, the levels of cross-polarized NRCS are larger than the co-polarized NRCS in the plane of ϕ s = 45 o , and the VV-polarized NRCS is the minimum.However, in the plane of ϕ s = 135 o , the behavior is the opposite.The behaviors of the aforementioned curves show a very good trend with conclusions drawn in other relevant papers [47,48].
Remote Sens. 2019, 11, x FOR PEER REVIEW 13 of 18 specular direction for co-polarized cases, which is a logical result.However, for cross-polarized results, the maximum value does not appear in the specular direction.In addition, unlike the backscattering configuration, the NRCS results for vertical-transmit-horizontal-receive (HV) polarization are not exactly the same as the horizontal-transmit-vertical-receive (VH) polarized results except for the specular direction.Furthermore, the NRCS results in fully bistatic configurations at Ka-band are investigated as well.Figure 11b,c  , the behavior is the opposite.The behaviors of the aforementioned curves show a very good trend with conclusions drawn in other relevant papers [47,48].

Normalized Radar Cross Section Dependence on Wind Speed
When referenced to 2D surfaces, it is of great interest to discuss the dependence of NRCS on the wind speed.In Figure 12, NRCS results which vary with wind speed are shown under forwardbackward configuration at Ka-band.From Figure 12, several items of importance may be deduced.First, the maximum of energy for co-polarizations is received around the specular direction, which is a logic statement.Second, the magnitude of NRCS for cross-polarizations increases with the increase

Normalized Radar Cross Section Dependence on Wind Speed
When referenced to 2D surfaces, it is of great interest to discuss the dependence of NRCS on the wind speed.In Figure 12, NRCS results which vary with wind speed are shown under forward-backward configuration at Ka-band.From Figure 12, several items of importance may be deduced.First, the maximum of energy for co-polarizations is received around the specular direction, which is a logic statement.Second, the magnitude of NRCS for cross-polarizations increases with the increase in wind speed.However, when it comes to the co-polarization modes, the scattering coefficient decreases with wind speed around the specular direction, whereas it increases with wind speed when the scattering direction is away from the specular direction.

Normalized Radar Cross Section Dependence on Wind Speed
When referenced to 2D surfaces, it is of great interest to discuss the dependence of NRCS on the wind speed.In Figure 12, NRCS results which vary with wind speed are shown under forwardbackward configuration at Ka-band.From Figure 12, several items of importance may be deduced.First, the maximum of energy for co-polarizations is received around the specular direction, which is a logic statement.Second, the magnitude of NRCS for cross-polarizations increases with the increase in wind speed.However, when it comes to the co-polarization modes, the scattering coefficient decreases with wind speed around the specular direction, whereas it increases with wind speed when the scattering direction is away from the specular direction.

Doppler Spectrum of Time-Evolving Sea Surfaces
In this part, Doppler spectrum characteristics at Ka-band is analyzed with the standard spectral estimation technique.In the following simulations, the normalized Doppler spectrum are averaged from 100 samples of sea surfaces with the size of 50 m × 50 m, and the wind speed of 5 m/s.The time step is set as 0.002 s and the total evolution time is 2.048 s.
Figure 13 compares the Doppler spectrum for up-wind direction with different incident angles.As indicated in Figure 12, it is observed that at moderate incident angles, the Doppler spectrum becomes wider compared with a lower incident angle.However, the variation of Doppler central frequency is not very obvious.In addition, the Doppler spectrum at Ka-band for different wind directions and polarizations are shown in Figure 14.In fact, researches have reported that there are two peaks in the cross-wind direction [5,7,31], which can be attributed to the fact that both the waves moving toward and away from the radar can be detected.However, the two peaks aliased to one peak due to the broader shape of the spectrum for Ka-band.

Doppler Spectrum of Time-Evolving Sea Surfaces
In this part, Doppler spectrum characteristics at Ka-band is analyzed with the standard spectral estimation technique.In the following simulations, the normalized Doppler spectrum are averaged from 100 samples of sea surfaces with the size of 50 m 50 m  , and the wind speed of 5 m/s.The time step is set as 0.002 s and the total evolution time is 2.048 s. Figure 13 compares the Doppler spectrum for up-wind direction with different incident angles.As indicated in Figure 12, it is observed that at moderate incident angles, the Doppler spectrum becomes wider compared with a lower incident angle.However, the variation of Doppler central frequency is not very obvious.In addition, the Doppler spectrum at Ka-band for different wind directions and polarizations are shown in Figure 14.In fact, researches have reported that there are two peaks in the cross-wind direction [5,7,31], which can be attributed to the fact that both the waves moving toward and away from the radar can be detected.However, the two peaks aliased to one peak due to the broader shape of the spectrum for Ka-band.
Figure 15 illustrates the comparison of Doppler spectrum for different polarizations.When the incident angle is 40 , it is observed that the Doppler shift at HH polarization mode is larger than that at VV polarization mode, and the Doppler shifts at cross polarization modes are the smallest.As for the Doppler width, the discrepancy caused by polarizations are relatively small.For the incident angle of 60 , the difference of Doppler width between HH and VV polarizations are obvious.However, there is almost no difference for the Doppler spectrum between VV and cross polarizations.That is to say for large incident angles, both the Doppler central frequency and the width are close for VV and HV polarizations.   Figure 15 illustrates the comparison of Doppler spectrum for different polarizations.When the incident angle is 40 • , it is observed that the Doppler shift at HH polarization mode is larger than that at VV polarization mode, and the Doppler shifts at cross polarization modes are the smallest.As for the Doppler width, the discrepancy caused by polarizations are relatively small.For the incident angle of 60 • , the difference of Doppler width between HH and VV polarizations are obvious.However, there is almost no difference for the Doppler spectrum between VV and cross polarizations.That is to say for large incident angles, both the Doppler central frequency and the width are close for VV and HV polarizations.

Conclusions
In this paper, FPFSM shows quite a good accuracy and computing efficiency in evaluating the radar returns from the electrically large sea surface, and it works well in a wide range of frequencies.In addition, the scattering characteristics from sea surface are predicted by FPFSM at Ka-band in bistatic configurations as a function of scattering angle, polarization mode, and wind speed.One can draw the following conclusions: (i) The maximum of energy is received around the specular direction in the forward-backward configuration for co-polarizations, but it does not hold true for cross polarizations.(ii) In cross polarization cases, the level of NRCS results increase with wind speed in both specular and diffuse scattering directions, which is different in the co-polarized cases.(iii) As for the Doppler spectrum, the results for different wave bands show similar characteristics in upwind and co-polarizations cases.However, for cross-wind direction, the Doppler spectrum has two peaks at low wave bands but the two peaks aliased to one at Ka-band.All of the numerical results show the good performance of the proposed method, and it can be further applied to analyze the scattering from targets over sea surfaces and to the SAR imaging simulators, which will benefit the recognition and detection of targets in a maritime environment.

Figure 2 .Figure 2 .
Figure 2. Comparison of NRCS results between the proposed method and the original SSA-II method under bistatic configuration in the near specular regime.

Figure 2 .
Figure 2. Comparison of NRCS results between the proposed method and the original SSA-II method under bistatic configuration in the near specular regime.

Figure 3 .
Figure 3.The global and local coordinates on a facet.

Figure 3 .
Figure 3.The global and local coordinates on a facet.

Figure 4 .
Figure 4.The NRCS results calculated by full-polarized facet based scattering model (FPFSM) with different samples for u 10 = 5 m/s and the Elfouhaily spectrum as a function of the wavenumber.(a) NRCS results at HH polarization mode; (b) NRCS results at VV polarization mode; (c) NRCS results at HV and horizontal-transmit-vertical-receive (VH) polarization modes; (d) Elfouhaily spectrum.

Figure 5 .
Figure 5.Comparison of the NRCS results calculated by FPFSM and SSA-II at ultra high frequency (UHF) band.(a) Monostatic configuration; (b) Forward-backward configuration.

Figure 5 .
Figure 5.Comparison of the NRCS results calculated by FPFSM and SSA-II at ultra high frequency (UHF) band.(a) Monostatic configuration; (b) Forward-backward configuration.

Figure 5 .
Figure 5.Comparison of the NRCS results calculated by FPFSM and SSA-II at ultra high frequency (UHF) band.(a) Monostatic configuration; (b) Forward-backward configuration.

Figure 8 .
Figure 8.Comparison of the NRCS results calculated by FPFSM and SSA-II at L-band with different wind speeds.(a) The wind speed is 8 m/s; (b) The wind speed is 10 m/s.

Figure 7 . 18 Figure 7 .
Figure 7.Comparison of the NRCS results calculated by FPFSM, SSA-II and real SAR data at C-band.(a) Monostatic configuration; (b) Forward-backward configuration.

Figure 8 .
Figure 8.Comparison of the NRCS results calculated by FPFSM and SSA-II at L-band with different wind speeds.(a) The wind speed is 8 m/s; (b) The wind speed is 10 m/s.

Figure 8 .
Figure 8.Comparison of the NRCS results calculated by FPFSM and SSA-II at L-band with different wind speeds.(a) The wind speed is 8 m/s; (b) The wind speed is 10 m/s.

18 Figure 9 .
Figure 9.Comparison of Doppler spectrum results calculated by FPFSM and SSA-II at L-band for different incident angles and polarizations.

Figure 9 .
Figure 9.Comparison of Doppler spectrum results calculated by FPFSM and SSA-II at L-band for different incident angles and polarizations.

Figure 10 .
Figure 10.Performance of the FPFSM at Ka-band.

Figure 10 .
Figure 10.Performance of the FPFSM at Ka-band.
gives the NRCS with the incident direction interesting to find that the discrepancy between the VH and HV polarizations became wider in these configurations.In addition, the levels of cross-polarized NRCS are larger than the co-polarized NRCS in the plane of -polarized NRCS is the minimum.However, in the plane of

Figure 12 .
Figure 12.Variation of Bistatic NRCS with wind speeds evaluated by FPFSM for different polarization states.(a) HH and HV polarizations.(b) VV and VH polarizations.

Figure 12 .
Figure 12.Variation of Bistatic NRCS with wind speeds evaluated by FPFSM for different polarization states.(a) HH and HV polarizations.(b) VV and VH polarizations.

Figure 13 .
Figure 13.Comparison of Doppler spectrum at Ka-band for different incident angles and polarizations.(a) HH polarization; (b) VV polarization; (c) HV polarization.

Figure 13 .
Figure 13.Comparison of Doppler spectrum at Ka-band for different incident angles and polarizations.(a) HH polarization; (b) VV polarization; (c) HV polarization.

Figure 15 .
Figure 15.Comparison of Doppler spectrum at Ka-band for different polarizations and incident angles.(a) Incident angle 40 ; (b) Incident angle 60 .

Table 1 .
Comparison of the time consumption (unit: seconds) between FPFSM and SSA-II.

Table 2 .
Comparison of the peak memory consumption (unit: Mb) between FPFSM and SSA-II.

Table 3 .
Mean errors and maximum errors (unit: dB) between FPFSM and SSA-II.

Table 1 .
Comparison of the time consumption (unit: seconds) between FPFSM and SSA-II.

Table 2 .
Comparison of the peak memory consumption (unit: Mb) between FPFSM and SSA-II.

Table 3 .
Mean errors and maximum errors (unit: dB) between FPFSM and SSA-II.