An FFT-Based Implementation of the First-Order Small-Slope Approximation for Simulating Microwave Scattering From Sea Surface

This article presents a simple and efficient implementation of the first-order small-slope approximation (SSA-1) to simulate the normalized radar scattering cross section (NRCS) from ocean surface based on fast Fourier transform (FFT) of the equivalent sea surface autocorrelation function. The implementation process only takes a pair of FFT operations between sea wave spectrum and sea surface autocorrelation function. Compared with the reported implementations of SSA-1 model based on the series expansion, the proposed method (called as “FFT-based implementation”) reduces the computation amount of FFT operations greatly. To examine the accuracy and efficiency of the FFT-based implementation relative to the direct integration implementation of SSA-1 (called as “classical implementation”), the monostatic and bistatic NRCSs simulated by using the two implementations with the help of Elfouhaily sea spectrum at various radar frequencies and sea states are compared. The comparison results indicate the FFT-based implementation not only holds almost the same accuracy as the classical implementation but also has a higher efficiency, especially for the simulation scenarios with small Rayleigh parameters. Furthermore, we discuss the influence of sampling range and sampling interval of sea spectrum components on the accuracy and efficiency of the FFT-based implementation and its feasibility on the simplified second-order SSA (SSA-2) model with the mathematic expression similar to SSA-1. The simulation results show that the sampling range and sampling interval of sea spectrum have different effect on the accuracy and efficiency of the FFT-based implementation, and the proposed implementation method is also feasible for the simplified SSA-2 model.

due to the impact of both local wind field and the long wave propagating from far away [2], which brings in many difficulties in modeling the radar scattering from ocean surface. Through decades of efforts by researchers, several analytical unified scattering models were developed for simulating the microwave ocean scattering, such as two scale model (TSM) [3], small-slope approximation (SSA) [4], integral equation model (IEM) [5], and its advanced version (AIEM) [6], and so on. Of these existing models, the SSA model seems more manageable mainly because, in the process of modeling ocean scattering, it uses the full wave components and also does not have to divide small-and large-scale roughness of sea surface by the artificial cut-off wavenumber as does TSM [3], [7]. The SSA model has been verified in a number of numerical simulations and also successfully used for interpretation of experimental data, e.g., [7], [8], [9], [10], [11], [12], [13], [14]. To date, the first-order SSA (SSA-1) and second-order SSA (SSA-2) models are usually appeared in literatures. However, considering the fact that the SSA-2 is only a small correction of its low-order SSA-1 but has a more cumbersome expression [7], [9], [11], the SSA-1 model is still widely used for the prediction of copolarized NRCSs from sea surface due to its relatively highcomputation efficiency and acceptable accuracy [13], [14], [15], [16], [17]. Therefore, in this article, we mainly focus on the implementation of SSA-1 model for simulating ocean microwave scattering.
The computational efficiency of SSA-1 model depends on its specific implementation methods. Typically, the SSA-1 model is implemented by using the direct numerical integration based on the discrete numerical solutions of the autocorrelation function of sea surface at different lag distance, e.g., [4], [7], [11], [14], [18]. Hereafter, the integration implementation is called as "classical implementation." In the process of the classical implementation, the numerical solutions of the sea surface autocorrelation function are obtained by integrating the sea wave spectrum over full wave components, and this step takes up most of the computing time for one NRCS simulation by SSA-1 with the classical implementation [14]. Hence, the computational efficiency of SSA-1 model based on the classical implementation is determined by the number of discrete samples of sea surface autocorrelation function required to be calculated to ensure computational convergence. The convergence is related to the characterization of integral kernel of SSA-1 model, which will be demonstrated in detail in Section II-B. For example, for low probing frequency and low wind speed (or say relatively small Rayleigh parameter), the classical implementation needs a larger integral domain and a finer discretization over it to ensure the accuracy of NRCS simulations [19], [20].
To improve the efficiency of SSA-1 in the simulation of NRCS at L-band and low wind speeds less than 3 m/s (i.e., at small Rayleigh parameter), the researchers expand the exponential term in SSA-1 expression as the Taylor series including higher order sea surface spectrum for the study of the transition from weak to strong diffuse specular scattering from sea surface [21]. The higher order sea surface spectrum contained in these expanded series terms are obtained from the sea surface autocorrelation function with different powers calculated by conducting FFT operations. This strategy improves the computation efficiency of SSA-1 model by using the same set of higher order sea surface spectrum evaluated at a given wind speed to calculate the NRCSs at different incidence angles. For the scenarios with small Rayleigh parameters in [21], a number of FFT operations are required to calculate the different higher order sea surface spectrum, e.g., about 150 times suggested and even more for the cases with larger Rayleigh parameters. To reduce the amount of FFT operations in the NRCS simulation, a relatively efficient implementation of SSA-1 by using only these effective terms of Taylor series with significant contribution to NRCS was proposed in our previous work [22]. Compared with the implementation based on series expansion [21], although using these effective terms defined in [22] can decrease the number of FFT operations and, thus, improve the computation efficiency to some extent, it still needs dozens of FFTs and even hundreds or thousands for the cases of large Rayleigh parameters, such as at Ku-band and high wind speed, which is not the desirable efficiency.
To further reduce the number of FFT operations in the implementations of SSA-1 model based on the Taylor series expansion mentioned above, this article proposes a faster implementation of SSA-1 by using an equivalent operation instead of using the Taylor series expansion. This work is the continuation of our previous work reported in [22]. The implementation proposed in this article, which hereafter is called as "FFT-based implementation," only needs two FFT operations, i.e., a pair of IFFT/FFT operations between the equivalent sea wave spectrum and the equivalent sea surface correlation function. Clearly, the FFTbased implementation can greatly reduce the number of FFT operations compared with the two existing implementations based on the series expansion [21], [22]. To examine the accuracy and efficiency of the FFT-based implementation, a comparison between the monostatic and bistatic NRCSs simulated by the FFT-based implementation and the corresponding simulations based on the classical implementation is conducted, aiming to extend an effective implementation method of SSA-1 model. The rest of this section is the organization of this article. In the following section, we introduce the methodology including microwave scattering geometry, the sea spectrum model, and the classical and FFT-based implementations of SSA-1 model. Next, the accuracy and efficiency of NRCS simulations by using the two implementations are compared in Section III. In Section IV, we further discuss the influence factors on the efficiency of the proposed FFT-based implementation as well as the feasibility of this implementation method on the simplified SSA-2 model. Finally, Section V concludes this article.

II. METHODOLOGY
A plane electromagnetic (EM) wave k i (θ i , φ i ) impinges upon a rough sea surface along the incidence angle θ i and azimuth angle φ i , and then a reflected wave k s (θ s , φ s ) with the scattering angle θ s and azimuth angle φ s can be observed. The radar measured scattering amplitude is associated with the sea surface roughness illuminated by antenna main beam, and the radar parameters, such as incidence angle, scattering angles, radar frequency, and polarization [1]. The general geometry of EM wave scattering from rough sea surface is illustrated in Fig. 1. Here, we define the x-axis and y-axis of the coordinate system point to upwind direction, i.e., φ w = 0 o , and crosswind direction, i.e., φ w = 90 o , respectively. In this coordinate frame, the negative direction of x-axis is in line with the downwind direction, i.e., φ w = 180 o .

A. Sea Spectrum Model
In scattering modeling, the roughness of the underlying sea surface shown in Fig. 1 is a key input for estimating the EM scattering from sea surface. Under the assumption that rough sea surface heights belong to the Gaussian distribution, the geometric description of sea surface can be expressed in terms of its power-density spectrum. In this study, the Elfouhaily power spectrum [23] is adopted due to its reliability in agreements with the slope model proposed by Cox and Munk [24] and the actual remote sensing data [7]. Generally, a two-dimensional (2-D) directional sea spectrum is expressed as the product of the omnidirectional spectrum S(K), representing the isotropic part, and the angular spreading function Φ(K, φ) [23]. The 2-D Elfouhaily spectrum is defined in polar coordinates as In (1) and (2), φ is the azimuth angle relative to upwind direction in wavenumber domain. K represents the wavenumber of sea wave component. Δ(K) is the upwind-crosswind ratio describing the anisotropy of the spreading function and given in [23]. The variance and autocorrelation function of sea surface elevation are directly related to the sea wave components by the following definitions: where ϕ is the azimuth angle relative to upwind direction in spatial domain.

B. Implementations of SSA-1 Model
The common form excluded the coherent components of SSA-1 is expressed as [4], [7], [11] where the subscripts p and q denote the received and transmitted polarization states, respectively; {k i , k s } are the horizontal projections of the incident and scattered wave vectors and {q i , q s } are the absolute values of their vertical projections, respectively. Q = q i + q s , and B pq (k s , k i ) is the polarization-related matrix of SPM-1 model and given in [7]. σ z is the root mean square height of rough sea surface. ρ(r) is the 2-D autocorrelation function of rough sea surface and r is the vector distance in spatial domain. For more detailed mathematical derivation of the SSA-1 model, readers are referred to [4], [7], [11]. 1) Classical Implementation: The classical implementation of SSA-1 model is based on numerical integration [4], [7], [11], [14]. In what follows, we briefly introduce the classical implementation reported in [11], where the integration variables in spatial domain are converted from 2-D vector r to the 1-D scalar r after completing the integration over azimuth directions. This conversion process is guided by rigorous mathematical derivation and finally is expressed as a closed form with the help of Bessel functions [11] where K B = |k s − k i |; J m (·) denotes the Bessel function of the first kind and of order m; I m (·) denotes the modified Bessel function of the first kind and of order m. φ w is the incident azimuthal angle relative to the upwind direction. ρ 0 (r) and ρ 2 (r) are the two parts of 2-D correlation function of sea surface ρ(r) in (4) with For integral implementation, one can see that the key is to obtain the discrete numerical solutions of correlation function by using (8), i.e., the integration over full sea wave spectrum components. This operation occupies most of the computation time for one NRCS simulation. The number of discrete samplings of correlation function required to be calculated depends on the behaviors of integral kernel in (6) shown at the bottom of this page. Taking the computation of the backscattered NRCS at L-band as an example, the variations of integral kernel function versus radar parameters and wind speeds are shown in Fig. 2.
It can be seen in Fig. 2 that it needs a longer integral distance and finer discretization to obtain an accurate NRCS simulation due to the strong oscillation of integral kernel function under the case of wind speed of 3 m/s and incidence angle of 60°, at which it corresponds to a smaller Rayleigh parameter. The discrete sampling numbers of autocorrelation function required to be calculated, i.e., ρ 0 (r) and ρ 2 (r) in (8), depend on the convergence and oscillation of integral kernel. In summary, the smaller the Rayleigh parameter is, the more computation time it requires before the integral kernel achieves to convergence.
2) FFT-Based Implementation: Let us go back to (5), the integral can be treated as the Fourier transform of an equivalent correlation function, C(r), i.e., the mathematic expression in curly braces in (5). Thus, (5) can be transformed into the following expression: where ρ(r) is the discrete 2-D autocorrelation function of sea surface expressed in (4) or (7). ρ(r) can be obtained by using the inverse FFT (IFFT) of the directional sea spectrum expressed in (1), which is the first FFT operation of the FFT-based implementation. Then, substituting the calculated ρ(r) into (10), a discrete 2-D equivalent correlation function C(r) is obtained. According to the Wiener-Khinchin theorem, the equivalent sea spectrum S E (K) can also be calculated from the equivalent correlation function C(r) by FFT operation, which is the second FFT operation. Thus, (9) can be further transformed to as follows: (12) Here, S E (k i − k s ) has the following form in rectangular and polar coordinates, respectively where k is the EM wavenumber; K B , which is determined by radar parameters and ranges from 0 to 2k, is a specific surface spectrum component along the azimuth of φ relative to upwind direction. According to abovementioned operations, for given radar parameters (i.e., radar frequency and observation geometry), the NRCS is directly related to the value of equivalent sea spectrum evaluated at the specific wavenumber of K B along the given wind direction of φ w relative to upwind direction. Hence, the critical work is to find the value of equivalent sea spectrum S E evaluated at the position of (K B , φ w ) in polar coordinates. Once the radar parameters and wind vector are given, the position of (K B , φ w ) is known. Since the S E calculated by FFT are organized in rectangular coordinate, the discrete 2-D matrix (K x , K y ) of spectrum variable in rectangular coordinates needs to be converted into its representations in polar coordinates. Thus, we get two corresponding 2-D matrixes of variables K and φ in polar coordinates. Next, we need search for these elements in matrixes K andφ whose values are close to K B andφ w , respectively. That is, the desired elements in matrixesK and φ need to meet the following conditions, respectively, where min { · } represents the operation of obtaining a minimum element in a matrix; Δ K and Δ φ are the threshold values, which both start from 0 and increase with the step of 0.1 until the elements satisfying (15) and (16) have the same indexes in two matrixes of K and φ, respectively. The same indexes, if exist, can be determined by the intersection operation between the two index arrays obtained from satisfying (15) and (16). If the same indexes have more than one, a compromise solution is adopted to select an unique index at which the corresponding element minimizes the summation of It should be note that it is easy for abovementioned operations to be conducted by using the matrix manipulation. Thus, an optimal value of S E (index) located at the specific K B and φ w is obtained to participate the subsequent NRCS calculation in the FFT-based implementation. The main work flows of abovementioned operations are shown in Fig. 3.

III. NUMERICAL RESULTS
In this section, we show the monostatic and bistatic NRCSs simulated at different radar frequencies, observation geometry, and sea states by using the FFT-based and classical implementations. Based on these NRCS simulations, the accuracy and efficiency of the two implementation methods are compared. In scattering modeling, the dielectric constant of sea water is function of radar frequency, seawater salinity, and sea surface temperature [25]. Here, we assume the sea surface temperature and salinity are constants of 20°C and 35‰, respectively. Besides, the fully developed sea with inverse wave age of 0.84 is assumed. For the classical implementation, to ensure the accuracy of simulation results as far as possible especially at small Rayleigh parameters and also to reduce variables, the sample interval of integral distance is uniformly set to be λ/200 for all the following NRCS simulations at different scenarios, where λ is the EM wavelength, and the required integral domain are determined by the convergence criterion that the absolute difference between the values of integral kernel calculated at several adjacent lag distances is less than 10 −8 . While, for the FFT-based implementation, the minimum sampling range of sea spectrum is set from 0 to 2k, where k is the EM wavenumber. The determination of this minimum sampling range is because  (15) and (16) in K and φ matrixes, respectively; the red dot is the determined index at which the value of the component of rough sea surfaceK B determined by radar parameters also varies between 0 to 2k for different radar observation configuration. The sampling interval of sea spectrum is ΔK = 2π/L, where L is the dimension of sea surface. In this study, to ensure the sampling accuracy of sea spectrum as much as possible, we set L = 200m, i.e.,ΔK = 0.03rad · m −1 . Therefore, the required sampling number N can be calculated by using the relationN = 2 × Round((2k · L)/(2π)), where the first factor of 2 results from the guarantee that the sampling sea spectrum components along all the azimuthal directions are between 0 to 2k. Based on above set of parameters, all the simulations are calculated by central processing unit (CPU) under the configuration: Intel(R) Core(TM) i9-10885H CPU @ 2.40 GHz. We can see in Fig. 4 that the NRCSs simulated based on the FFT-based implementation have a good agreements with that simulated with the classical implementation for both HH-and VV-polarizations. The root mean square error (RMSE) between the NRCSs simulated based on the two implementations at various wind speeds and radar frequencies are less than 0.3 dB. This means the FFT-based implementation well holds the simulation precision of SSA-1 model. Their computation efficiency is shown in Fig. 5.

A. Monostatic NRCS Simulations
At L-band, one can see from Fig. 5(a1) that the computation time for an NRCS simulation with the FFT-based implementation is about 6 s and it is independence on incidence angles and wind speeds. This is because the computation efficiency of FFTbased implementation is determined by the sampling number N , which is related to both the sampling interval, i.e., dimension of sea surface L, and the sampling range of sea spectrum. While, for the classical implementation, see Fig. 5(a2), the computation time is a few hundred seconds, which is dozens of times larger than that of the FFT-based implementation. Moreover, the computation time of classical implementation increases with the decrease of wind speed especially at relative large incidence angle. This is because the relatively small Rayleigh parameters due to the larger incidence angle and lower wind speed result in a longer integral domain of integral kernel (see Fig. 2).
At C-band, for the case of relatively high radar frequency shown in Fig. 5(b1), the computation time of the FFT-based implementation is about 120 s, which is larger than that of L-band. The reason is that the EM wavenumber k of C-band is larger than that of L-band, and thus, the number of samples within 0 to 2k with the same sampling interval increase. On the effect of the sampling range and sampling interval of sea spectrum on accuracy and efficiency of the FFT-based implementation, we will discuss it in Section IV. On the contrary, for the classical implementation, the computation time at C-band [see Fig. 5(b2)] is less than that of L-band [see Fig. 5(a2)] for different wind speeds, which results from the shorter integral domain of integral kernel due to the relatively large Rayleigh parameters at C-band. Besides, we can see from Fig. 5(c1) and (c2) that, compared with the classical implementation, the computation time for completing an NRCS simulation by the FFT-based implementation is greatly reduced, especially for L-band [see Fig. 5(c1)]. For example, at the scenario for L-band and low wind speed of 5 m/s, the efficiency has been improved by tens of times generally, and even more than a hundredfold increase at large incidence angles greater than 40°.
In addition to investigating the variation of NRCS versus incidence angles, we also examine the variation of NRCSs simulated by the two implementation methods with wind directions.  shows the corresponding simulation results at wind speed of 10 m/s and incidence angle of 40°for L-and C-bands. In this figure, we can see that the differences between the simulated NRCSs at different wind directions are all very small with the RMSE less than 0.1dB for both L-and C-bands. It suggests the FFT-based implementation is competent to simulate the dependence of NRCS on wind directions. The computation time of each one NRCS simulation is similar to the case of incidence angle shown in Fig. 5.

B. Bistatic NRCS Simulations
In above Section III-A, the accuracy of the FFT-based implementation in simulating the backscattered NRCSs at various radar frequencies, incidence angles, wind speeds, and wind directions are examined. To further validate the accuracy of this implementation for other radar observation geometries, we compare the bistatic NRCSs simulated at several general observation scenarios, which are shown in Fig. 7. Fig. 7(a1) and (a2) shows the specular NRCSs simulated by using the two implementation methods at L-and C-bands when wind speed is 10 m/s and wind direction is along the upwind direction. The simulated NRCSs versus scattering angles in incident plane are shown in Fig. 7(b1) and (b2) for L-and C-bands, respectively, when the incidence angle is 40°and wind speed is 10 m/s. The negative scattering angles in Fig. 7(b) indicate the scattering in the back region; while the positive ones represent the forward scattering. Besides, the NRCSs varying with scattering azimuth angles at wind speed of 10 m/s are also simulated at the same incidence and scattering angles of 40°A uthorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply. Fig. 7. Comparisons of simulated bistatic NRCSs by using the FFT-based and classical implementations at L-and C-bands along upwind direction for the cases of (a1) and (a2) scattering in specular direction, (b1) and (b2) scattering versus scattering angles in incident plane with incidence angle of 40 degrees, and (c1) and (c2) scattering versus scattering azimuth angles with the same incidence and scattering angles of 40°.
when the incident azimuth angle is along the upwind direction. The corresponding simulation results at L-and C-bands are shown in Fig. 7(c1) and (c2), respectively. The presence of the trough for HH-polarization is caused by the mathematical criteria used in the theoretical development of SSA-1 model, which has no real physical meaning [14], [18]. By comparisons, we can see these simulated bistatic NRCSs based on the two implementation methods have the very small difference with the RMSE less than 0.2 dB. The good consistency indicates the FFT-based implementation used for simulating bistatic scattering is also feasible.

IV. DISCUSSION
In Section III, we mention that the efficiency of the FFT-based implementation is only related to the number of samples participating in FFT operation. However, the number of samples are determined by the sampling interval of ΔK, i.e., the dimension of sea surface L, and the range of sea spectrum components from 0 rad/m to its upper limit K max , i.e., (0, K max ). Therefore, in this section, we further investigate the effect of L and K max on accuracy and efficiency of the FFT-based implementation as well as the feasibility of this implementation method on the simplified SSA-2 model.

A. Influence Factors on Accuracy and Efficiency of the FFT-Based Implementation
The simulated NRCSs in Section III are conducted by the FFT-based implementation with the parameter combination of K max and L equal to(K max = 2k, L = 200 m). Hereafter, we use the combination of(K max , L) to represent the values of corresponding parameters of sampling range and sampling interval of sea spectrum used in the FFT-based implementation. Therefore, to investigate the effect of K max and L on the FFT-based implementation, it needs to compare the differences between the simulated NRCSs under the different combinations of (K max , L). In what follows, we take the backscattered NRCSs simulated by using the FFT-based implementation with the combination of (2k, 200) in Fig. 4 as the reference for comparisons. For comparison, we need to simulate the backscattered NRCSs by using the FFT-based implementation with other two combinations of (2k, 100) and (3k, 200) at wind speeds of 5-, 10-, and 15-m/s. The corresponding comparisons of VV-polarized NRCS at L-and C-bands are shown in Figs. 8 and 9, respectively. The HH-polarization results are not shown because of its similar behaviors to that of VV-polarization. Besides, the absolute errors between them for an NRCS simulation at each incidence angle are also presented at the right-hand ordinate axis (marked in red) in these figures.
We can see from the left columns in Figs. 8 and 9 that the simulated NRCSs at L-and C-bands by using the FFT-based implementation at various incidence angles for the three wind speeds have almost no difference between the combinations of (2k, 200) and (3k, 200), which means these sea wave components with the wavenumber larger than 2k has almost no effect on NRCS simulations. While, as shown in the right columns of Figs. 8 and 9, there are relatively obvious difference between the results simulated under the combinations of (2k, 200) and (2k, 100), and the difference becomes large with the increasing wind speed. The comparisons in Figs. 8 and 9 indicate that the parameter of L has a relatively large influence on NRCS simulations, that is, the FFT-based implementation is more sensitive to the sampling interval of sea spectrum. A smaller value of L can lead to a larger sampling interval of sea spectrum, which may omit some low-frequency wave components with significant contribution to sea surface roughness. When these low-frequency wave components near the peak frequency of sea spectrum are omitted in the process of sea spectrum component sampling, the sea surface roughness, such as variance of sea surface height, can be affected and more significantly in high sea states. This is the reason why different sizes of L can cause the difference on the NRCS simulation. To check the influence of sampling interval on sea surface variance, we calculate the variances of sea surface elevation at the three different wind speeds by using the FFT method represented in (4) under the three different combinations of (K max , L). For comparison, we also calculate the variances of sea surface elevation by using the numerical integral method represented in (3) over full wave components. The corresponding comparisons are shown in Fig. 10.  In Fig. 10, we can see that these high-frequency wave components with the wavenumber ranging from 2k to 3k have almost no influence on variance of sea surface height for different wind speeds when the sampling interval is fixed; while the sampling interval of sea spectrum, i.e., L, has a different effect on the variance of sea surface height under different wind speeds. The effect of sampling interval is significant for high wind speed, e.g., 15 m/s, but a little for the low-and medium-wind speeds. This may be because the peak of sea wave spectrum moves toward the direction of a smaller wavenumber with the increase of wind speed [23]. Thus, a large sampling interval (or inadequate sampling) may lead to the absence of some large-scale wave components, in particular, the absence of wave components corresponding to the peak of the sea spectrum, in the reconstruction of sea surface roughness. These missing large-scale long wave components can degrade some accuracy, although limited amount, of the NRCS simulations, which can be seen from the comparison between the simulated NRCSs by using L = 100 m and L = 200 m. Unlike the effect of L and K max on sea surface roughness parameters, however, the two parameters both have a remarkable effect on computation time of FFT-based implementation method because they directly determine the number of samples used in FFT operation. The comparison of computation time based on the FFT-based implementation under the three combinations of (K max , L) at wind speed of 5 m/s for L-and C-bands are shown in Fig. 11, respectively. As is shown in Fig. 11 that, compared with the combination of (2k, 100), the computation time based on the FFT-based implementation increases when the values of K max or L become large.

B. Determination of the Minimum Size of Sea Surface Length
By the analysis in Section IV-A, we know that the size of sea surface length has influence on both the accuracy and efficiency of the simulation results, but more influence on the efficiency. Therefore, an appropriate sea surface size is very important to improve the efficiency of the FFT-based method. To preserve the true statistical characteristics of sea surface as much as possible, we need to sample the sea wave components, K p , at which the sea spectrum reaches to its maximum. Thus, it requires that K p should be an integral multiple of the sampling interval, ΔK = 2π/L. That is where n is a positive integer. Moreover, K p of Elfouhaily sea spectrum can be correlated with wind speed and expressed as [23] where Ω is the inverse wave age and set to be 0.84 in this article; g = 9.81 m/s 2 is the gravity acceleration; U 10 represents the wind speed at the height of 10 m above the sea surface.
Combining (17) and (18), the relationship between wind speed and sea surface length becomes neat Equation (19) ensures that the wave component with the highest energy can be sampled in the sampling process of sea spectrum, preserving the rough statistical characteristics of sea surface as much as possible, such as variance. The longer the sea surface length is, the smaller the sampling interval of sea spectrum is, resulting in an increase in the number of sampling points and the degradation in efficiency. Therefore, we need to input a sea surface length as small as possible in the algorithm to improve its efficiency.
Given this, we examine the statistical characteristics of rough sea surface and the accuracy of scattering simulation when n is set to be 1 (i.e., n = 1). Thus, L = 0.91 · U 2 10 is used to determine the length of sea surface adopted in the simulations of sea surface statistical characteristics and microwave scattering at different sea states. Similar to Figs. 10 and 12 shows the results of sea surface variance calculated by (3) with full wavenumber components and by (4) with the sampled wave components based on the combination of (2k, L = 0.91 · U 2 10 ) at L-band for different wind speeds.
It can be seen that the geometric statistical characteristics of sea surface can be preserved effectively by IFFT of the sampled sea spectrum based on the selection criterion of sea surface length, i.e., L = 0.91 · U 2 10 . Under this criterion, the number of sampling points for FFT at low wind speed is greatly reduced relative to the fixed set of sea surface length previously, L = 200 m or L = 100 m, resulting in a higher efficiency for the FFT-based implementation.
To further examine the effect of this criterion on the accuracy of the FFT-based implementation at different wind speeds and radar frequencies, we compared the backscattered NRCSs   Fig. 13.
One can see in Fig. 13 that differences of simulated NRCSs based on the two combinations of (2k, 200) and (2k, 0.91 · U 2 10 ) are small and the maximum RMSE is less than 0.1 dB. This indicated that the minimum sea surface length by using the criterion, L = 0.91 · U 2 10 , can also maintain the high accuracy for scattering simulation. It is worth noting that using L = 0.91 · U 2 10 can reduce the number of samples in the FFT operation and improve the efficiency of the FFT-based implantation especially for low and medium wind speeds, which can be seen clearly in Fig. 14. For example, at low wind speed of 5 m/s, the runtime for an NRCS simulation can be further reduced from about 6 s and 120 s under the combination of (2k, 200) to 1 s and 2.5 s for L-band and C-band, respectively. This is not a small efficiency gain.

C. Feasibility of the FFT-Based Implementation on SSA-2
According to the abovementioned analysis, the FFT-based implementation is efficient and appropriate for the scattering model with a mathematic expression similar to that of SSA-1 represented in (5). Therefore, if the SSA-2 can be simplified with the mathematic expression similar to (5), the FFT-based implementation should be applicable to SSA-2 in theory. Fortunately, to facilitate the calculation of SSA-2 model, Voronovich et al. [7] make a legitimate approximation for the small correction term in SSA-2, based on which the expression of SSA-2 model is degraded into the one similar to SSA-1 after the correlation function of sea surface is substituted by a modified correlation function. The final mathematic expression of SSA-2 model with the modified correlation function is expressed as [7] with Here, h M (r) is the modified correlation function of sea surface. M pq (k s , k i ; ξ) are related to polarization and their explicit expressions can be seen in [7]. The remaining parameters are the same as (5). The corresponding modified 2-D sea spectrum is where S(K) is 2-D sea spectrum expressed as in (1). For more details about the simplification process, readers can be referred to [7]. With the help of the modified sea spectrum in (22), the FFT-based implementation can be applied to the simplified SSA-2 model.
To examine the simulation accuracy of the SSA-2 model with the FFT-based implementation, the backscattered NRCSs simulated by using the simplified SSA2 based on the integral method at C-and Ku-bands along upwind direction were extracted from [7] to serve as the reference data for comparison. There are no the reference data for HH-polarization for C-band in [7]. For C-band, we used the minimum sea surface length determined by L = 0.91 · U 2 10 in the NRCS simulations. While, for Ku-band, considering the burden of computing memory resulting from the huge number of samplings of sea spectrum at the range of 0 to 2k, we selected a fixed value of L, i.e., combination of (2k, 100 m), to execute SSA-2 for simulating the NRCSs at Ku-band. The corresponding simulation results are shown in Fig. 15.
In this figure, an overall good consistency between the simulated NRCSs based on the FFT-based implementation and the reference data generated by the integral implementation is observed except for Ku-band at wind speed of 15 m/s. The difference may come from the imprecise sampling of sea spectrum due to the smaller L (here, L = 100 m). Overall, the FFT-based implementation is feasible for SSA-2 to some extent. The runtime for C-band simulations by using SSA-2 is about the same to that of SSA-1 model (see Fig. 14); while the runtime for an NRCS simulation at Ku-band by using SSA-2 is about 320 s under the combination of (2k, 100 m).
Considering the computational time that it requires for further improving the simulation accuracy of FFT-based implementation through adding more dense samplings at Ku-band, the classical implementation is perhaps easier to operate for the ocean scattering simulation at high probing frequencies. This is because the large Rayleigh parameters at high frequency can shorten the integral domain to be required [11], [18].
To sum up, compared with the classical implantation, the advantage of the FFT-based implementation is that it has a higher computational efficiency in the simulation of NRCSs from sea surface at low and medium radar frequencies, especially for L-band under the cases of low wind speeds (i.e., relatively small Rayleigh parameters).

V. CONCLUSION
This article presents an effective implementation of SSA-1 model based on two FFT operations: the first one is the IFFT from sea spectrum to autocorrelation function of rough sea surface, and the other is the FFT from the equivalent correlation function of sea surface to the equivalent sea spectrum. This implementation greatly reduces the amount of FFT operations compared with the existing implementations based on series expansion reported in [21] and [22] where they both need dozens or even hundreds of FFTs. To check the accuracy and efficiency of the proposed FFT-based implementation in the simulation of ocean microwave scattering, we take the simulation results by using the classical implementation as the reference for comparison. The mono-and bistatic NRCSs at L-and C-bands are first simulated based on the FFT-based implementation under the combination of (K max , L) equal to (2k, 200 m) and the classical implementation at different wind speeds of 5-, 10-, and 15-m/s. Then, the corresponding simulation results with the two implementations are compared to evaluate the accuracy and efficiency of the FFT-based implementation. Furthermore, the influence of the sampling range of sea spectrum, (0, K max ), and its sampling interval, 2π/L, on the accuracy and efficiency of the FFT-based implementation along with the feasibility of this implementation method on the simplified SSA-2 model are also discussed. Based on these NRCS simulations, some conclusions are summarized as follows.
1) The NRCSs simulated by SSA-1 model with the two implementations have a good consistency with the RMSE less than 0.3 dB, and the efficiency of the FFT-based implementation is relatively high compared with the classical implementation especially at low-and medium-frequency radar bands. 2) These sea wave components with wavenumber lager than 2k has very little influence on simulation accuracy, but the value of L has a relatively large effect on the accuracy of the FFT-based implementation for the case of high wind speeds. For the computational efficiency of FFT-based implementation, it is obviously influenced by both K max and L, especially for high probing frequencies, e.g., Ku-band. The advantage of the proposed method in this article exists in its efficiency in the simulation of ocean microwave scattering at low and medium probing frequencies by using the sea surface length determined by L = 0.91 · U 2 10 . For the case of high-frequency and high wind speed, the classical implementation method may be more manageable. Therefore, the FFT-based and the classical implementations of SSA-1 model can complement each other. It is advisable to choose an appropriate implementation method for the specific situations.
3) The FFT-based implementation can also be transplanted on the simplified SSA-2 model with the mathematic expression similar to that of SSA-1 within an acceptable accuracy. Overall, the contribution of this article is to provide another efficient implementation, which complements the classical implementation of SSA-1 model each other in the simulation of microwave scattering from sea surface.