# A Novel Hybrid Algorithm for Transient Near-Field Scattering from Rough Surface in 2-D Case

Authors: Wei Tian; Bing Wei; Qian Yang
Journal: International Journal of Antennas and Propagation
 (2021)

---

## Abstract

A novel hybrid algorithm is proposed to reduce the computation cost of the finite-difference time-domain (FDTD) method in calculating the transient near-field scattering from rough surface. The scattering problem is split into the FDTD calculation of equivalent sources on the contour enclosing rough surface and the calculation of the near-field radiation with the two-dimensional (2-D) time-domain Huygens’ principle. The radiation fields are found from a surface integral of the temporal convolution for which the direct numerical integration of the convolution is computationally expensive. In this paper, the 2-D time-domain Green’s function as the convolution kernel is approximated with a sum of exponential terms by using the Prony’s method. Then, the semianalytical recursive convolution (SARC) approach is applied to complete the update of the near-field radiation. Compared with the traditional FDTD, this hybrid algorithm can significantly reduce the memory usage and run time, especially for the large distance between the rough surface and observation point.


---

## Body

## 1. Introduction

With the extensive application of the wide-band radar, the transient electromagnetic (EM) scattering from the rough surface has gained increasing attention, which plays important role in the field of radar surveillance, EM imagining, target identification, and stealth, etc. [1–4]. The previous researches are mainly focused on the simulation of the EM scattering from rough surface in the far-field zone [5–7]. In fact, the near-field scattering problem is the widespread presence of the target working above background, such as the low-altitude cable over ground or sea. Only by mastering the near-field scattering characteristics from rough surface, we can systematically analyze the EM interference from background and obtain the relatively complete coupling result of target. Therefore, it is of critical importance to study an efficient time-domain algorithm on the near-field scattering for the coupling analysis of the low-altitude target and the protection of the sensitive electronic device.In the past decade, the numerical method has been widely used in the EM analysis owning to the rapid development of the high performance computer [8–11]. Among them, the finite‐difference time-domain (FDTD) method is a powerful and popular technique for numerically solving the transient field problem, and it has achieved successful applications in the far-field scattering from rough surface [12–14]. However, when the FDTD is used to simulate the near-field scattering, a brute-force way is to enlarge the computational domain to enclose the rough surface and observation point, which will lead a heavy burden to the computation resource available. Moreover, the farther the near-zone observation point is from rough surface, the heavier the computation burden will be. By now, some hybrid algorithms combining the FDTD and analytical-approximate method have been proposed to reduce the FDTD space lattice and improve the computational efficiency. In [15], the FDTD associated with the Kirchhoff surface integral was proposed to calculate the near-field response of the dipole antenna in the three-dimensional (3-D) case, where the dipole antenna was assigned in the FDTD region, and the near-field response outside the FDTD space was calculated by the Kirchhoff surface integral. For the two-dimensional (2-D) problem, by means of the 2-D time-domain Huygens’ principle, the temporal convolution integral needs to be solved to obtain the near-field response exterior to the FDTD region, which complicates the computation even more. The direct integration of the temporal convolution is computationally expensive, and hence it is applicable for the EM analysis of the small model such as the infinitely long line source in [16, 17]. Xu et al. in [18] employed the trapezoidal approximation to solve the temporal convolution, and then the near-field scattering exterior to the FDTD region was computed by applying the scheme of treating the radiation field as the equivalent cylindrical wave, which can be accurately used for the target with cylinder shape. However, for the large rough surface model, the characteristics of the structure and scattering wave make the simulation of the near-field scattering with these approaches difficult. In recent years, a fruitful semianalytical recursive convolution (SARC) approach was introduced in the dispersive medium problem, realizing the high-efficient computation of the convolution integral [19, 20]. But when the 2-D time-domain Green’s function is acted as the convolution kernel, the recursive operation is difficult to realize due to its specific structure. In this paper, this difficulty is overcome by applying the Prony’s method to approximate the convolution kernel with a sum of exponential terms, which avoids the complex analysis of the convergence of the series expansion scheme in [21]. The hybrid algorithm combining the FDTD with the 2-D time-domain Huygens’ principle is proposed to calculate the transient near-field scattering from the rough surface, in which the normal part of the convolution integral is calculated by the SARC approach and the singular part is done using the linear interpolation. Because a great deal of FDTD meshes between the rough surface and near-zone observation point are removed, the memory usage and run time in the proposed algorithm are dramatically reduced than those in the traditional FDTD method. Meanwhile, this hybrid algorithm retains the unique advantages of the FDTD method.The rest of this paper is organized as follows: in Section2, the hybrid algorithm combining the FDTD with the 2-D time-domain Huygens’ principle is elaborated, including the Prony approximation of the convolution kernel and the SARC approach for solving the near-field radiation of the equivalent sources. In Section 3, the accuracy and efficiency of the proposed hybrid algorithm are verified through the comparison of the numerical results by the hybrid algorithm and by the traditional FDTD method. Section 4 draws the conclusions of this paper and offers the recommendations for future work.


## 2. Hybrid Algorithm for Near-Field Scattering from Rough Surface

For the FDTD calculation of the transient near-field scattering, a common way is to expand the FDTD region to enclose the rough surface and the near-zone observation point, as shown in Figure1(a). With the increase in the height of observation point, the memory usage and run time will sharply increase, thus often making it difficult to realize the FDTD calculation of the near-field scattering. For the proposed hybrid algorithm, however, the FDTD region only contains the rough surface as indicated in Figure 1(b); the near-field response at the observation point P is regarded as the radiation of the equivalent sources on the output boundary of the FDTD region. Comparing Figures 1(a) and 1(b), it can be seen that the computational domain in the hybrid algorithm is much smaller than that in the FDTD method, especially in the case with high observation point.Figure 1
Comparison of the computational domain for the near-field scattering. (a) Traditional FDTD method. (b) Proposed hybrid algorithm.
(a)(b)


### 2.1. Theory of Near-Field Transformation in 2-D Case

Assume thatE and H are respectively the electric and magnetic fields on the output boundary surface of the FDTD domain. According to the equivalence principle, the surface electric and magnetic currents J,M and charges ρs,ρm can be obtained by(1a)J=n^×H,M=−n^×E,(1b)n^⋅μH=ρs,n^⋅μE=ρm,where n^ denotes the unit normal vector on the output boundary. For the 2-D problem, the vector and scalar potentials can be expressed with the 2-D time-domain Green’s function as(2a)Ar,t=μ4π∫Ldl′∫0t−R/cJr′,t′t−t′2−R2/c2dt′,(2b)Fr,t=ε4π∫Ldl′∫0t−R/cMr′,t′t−t′2−R2/c2dt′,(2c)φr,t=14πμ∫Ldl′∫0t−R/cρsr′,t′t−t′2−R2/c2dt′,(2d)φmr,t=14πμ∫Ldl′∫0t−R/cρmr′,t′t−t′2−R2/c2dt′.Here,μ and ε represent the permeability and permittivity, R=r−r′ is the distance from the observation point r to source field r′, and c is the speed of the EM wave. In accordance with the vector and scalar potentials, the 2-D radiation field exterior to the FDTD region can be written as(3a)Er,t=−1ε∇×Fr,t−∇φ−∂∂tAr,t,(3b)Hr,t=1μ∇×Ar,t−∇φm−∂Fr,t∂t.Substituting equation (2) into equation (3) and shifting the time reference backward by R/c, the formals for calculating the 2-D radiation field of the equivalent sources can be obtained from [21], in which, for the TM case, are(4a)Ezr,t+Rc=12π∫Ldl′−μ∫0t∂Jzr′,t′/∂t′t+R/c−t′2−R2/c2dt′−y−y′R2∫0tt+R/c−t′∂Mxr′,t′/∂t′t+R/c−t′2−R2/c2dt′+x−x′R2∫0tt+R/c−t′∂Myr′,t′/∂t′t+R/c−t′2−R2/c2dt′,(4b)Hxr,t+Rc=12π∫Ldl′−ε∫0t∂Mxr′,t′/∂t′t+R/c−t′2−R2/c2dt′−x−x′μR2∫0tt−t′∂ρmr′,t′/∂t′t+R/c−t′2−R2/c2dt′+y−y′R2∫0tt+R/c−t′∂Jzr′,t′/∂t′t+R/c−t′2−R2/c2dt′,(4c)Hyr,t+Rc=12π∫Ldl′−ε∫0t∂Jzr′,t′/∂t′t+R/c−t′2−R2/c2dt′−y−y′μR2∫0tt+R/c−t′∂ρmr′,t′/∂t′t+R/c−t′2−R2/c2dt′+x−x′R2∫0tt+R/c−t′∂Jzr′,t′/∂t′t+R/c−t′2−R2/c2dt′.The corresponding formulas in the TE mode can be derived with the principle of duality, which are not described here. It is obvious that the convolution kernel in equation (4) has a square-root singularity at t′=t. An efficient way in addressing the singularity is to approximately calculate the convolution integral in the interval t−Δt,t using the linear interpolation. Assuming t−t′=τ, the time derivative of the equivalent source, e.g., Mr,t′, in the time-reversed and offset form can be given by(5)MnΔt−τ≈Mn+Mn−1−MnΔtττ∈0,Δt.For instance, substituting (5) into the second term of (4a), we have(6)∫t−Δttt+R/c−t′Mr′,tt+R/c−t′2−R2/c2dt′=Mn+RMn−Mn−1cΔtτ+Rc2−R2c2|τ=0τ=Δt+Mn−1−MnΔtτ+R/c2τ+Rc2−R2c2+R22c2lnτ+Rc+τ+Rc2−R2c2τ=0τ=Δt.A similar procedure can be implemented for other terms in equation (4) to obtain the corresponding approximate results in the interval t−Δt,t. The convolution integral in the interval 0,t−Δt is computed by using the high-efficient SARC approach, as will be elaborated later.


### 2.2. Prony Approximation of Convolution Kernel

Although equation (4) is applicable for the computer execution, the direct integration of it is such cumbersome as to be impractical, due to the need for storing and processing all past values of the time derivative of the equivalent source. The computational complexity of this scheme scales as ONmax2, where Nmax is the total time step. In this paper, the SARC approach will be applied to decrease the complexity of the convolution integral. It can be found from equation (4) that the convolution kernel mainly consists of two forms: 1/t+R/c2−R2/c2 and t+R/c/t+R/c2−R2/c2, which make it difficult to execute the recursive operation of the convolution. Based on the idea of the signal estimation in the high-speed circuit [22], the Prony’s method is used to approximate the convolution kernel with a sum of the exponential terms:(7)fN≈∑i=1QCiexpAiN=∑i=1QCiμiN.Here,fN is the discrete sample of the convolution kernel at t=NΔt, where N denotes the equally spaced time step, Q represents the total number of poles, and Ai and Ci represent the coefficient and amplitude of the ith pole. The Ai and Ci are evaluated in the Prony’s method by solving two sets of linear equations with an intermediate nonlinear equation with Qth degree. Firstly, Nm equally spaced time samples of the convolution kernel are taken to set up a set of linear equation as(8)f0f1f2⋯fQ−1f1f2f3⋯fQ⋯⋯⋯⋱⋯fNm−Q−1fNm−QfNm−Q+1⋯fNm−2α1α2⋯αQ=−fQ−fQ+1⋯−fNm−1,where Nm−Q>Q. The roots α1,α2,…,αQ can be determined by solving this oversampled set of linear equations with the least-squares approach. By virtue of equation (8), the μ can be found as the roots of the following Qth degree equation:(9)μQ+αQμQ−1+αQ−1μQ−2+⋯+α2μ+α1=0.If the complex conjugate roots arise in the solution to equation (9), use their magnitude to substitute for the complex pair of μ, and cancel any duplicate. Substituting μi to equation (7) and introducing Nm discrete samples of the convolution kernel, a new set of the linear equations with the Q values of the Ci can be established as follows:(10)111⋯1μ1μ2μ3⋯μQμ12μ22μ32⋯μQ2⋯⋯⋯⋱⋯μ1Nm−1μ2Nm−1μ3Nm−1⋯μQNm−1C1C2C3⋯CQ=f0f1f2⋯fNm−1,where Nm>Q. The coefficients C1,C2,…,CQ can be obtained by using the least-squares approach to equation (10).Figure2 shows the comparison of the fitting results by the Prony’s approximation with the analytical results, in which the red short dash denotes the analytical result and blue line denotes the fitting result with Prony’s method. For the fitting calculation, the space increment is Δ=0.05m, the time increment is Δt=Δ/2c, and the observation point and source point are located at 20Δ,80Δ and 30Δ,10Δ, respectively. The sampling of the convolution kernel begins at intR/cΔt+1, and the interval is 20Δt. In Figures 2(a) and 2(b), the 1/τ+R/c2−R2/c2 and τ+R/c/τ+R/c2−R2/c2 are respectively expressed as convolution kernels 1 and 2 for convenience. The total numbers of the samples and poles are selected to be 30 and 10, respectively. It can be observed from Figure 2 that the nearly identical results are produced with the two methods, implying the reliability of the Prony’s method.Figure 2
Comparison of the fitting results by the Prony’s approximation with the analytical results. (a) Convolution kernel 1. (b) Convolution kernel 2.
(a)(b)


### 2.3. Semianalytical Recursive Convolution Approach

The RASC approach is a fruitful scheme for the calculation of the convolution integral, which is based on the idea of the efficient analysis of the transfer response in the digital signal processing (DSP) [23]. For a linear time-invariance and dynamical system, as shown in Figure 3, the system response yt can be expressed as the convolution of the input signal xt with the impulse response ht. We have(11)yt=ht∗xt=∫0tht−t′xt′dt′.Figure 3
Transfer response for linear time-invariance system.In general, the impulse response in an asymptotically stable system is of the form,(12)ht=∑q=1QHqexp−βqtUt,where q stands for the number of poles, and Ut is the switch function:(13)Ut=1t≥00t<0.Substituting equation (12) into equation (11), the system response at t=nΔt can be obtained as(14)yn=∑q=1Q∫0nΔtHqexp−βqnΔt−t′xt′dt′.The system response for theqth pole is(15)yqn=∫0nΔtHqexp−βqnΔt−t′xt′dt′.Dividing the integral of equation (15) into two intervals, 0,n−1Δt and n−1Δt,nΔt, we get to the recursive formula of system response as(16)yqn=exp−βqΔtyqn−1+Hqexp−βqnΔt∫n−1ΔtnΔtexpβqt′xt′dt′.In this paper, the near-field radiated by the equivalent source is analogous to the system response in the DSP. The convolution kernel can be regard as the impulse response after it is approximated with a sum of exponential terms. The time derivative of the equivalent source corresponds to the input signal. Thus, the near-field radiation can be updated by applying the SARC approach. Next, focus on the SARC approach for evaluating the radiation field in equation (4), where the representative convolutions are of the forms,(17a)It+Rc=∫0t−ΔtWr′,t′t+R/c−t′2−R2/c2dt′,(17b)Kt+R/c=∫0t−Δtt+R/c−t′Wr′,t′t+R/c−t′2−R2/c2dt′.For equation (17a) as example, using the Prony’s method mentioned above can approximate the convolution kernel as(18)1t+R/c2−R2/c2≈∑q=1QCqexpAqt.According to equation (16), the temporal convolution in the interval 0,n−1Δt can be expressed as(19)Iqn−1=expAqΔtIqn−2+CqexpAqnΔt∫n−2Δtn−1Δtexp−Aqt′Wt′dt′,where Wt can be approximated in accordance with the linear interpolation approach as follows:(20)Wt′≈Wn−1+Wn−1−Wn−2Δtt′−n−2Δt.Substituting equation (20) into equation (19) and analytically solving the integral over the interval n−2Δt,n−1Δt, we find the following field-update equation:(21)Iqn−1=expAqΔtIqn−2+B0,qWn−1+B1,qWn−2,where(22)B0,q=−CqexpAqΔtAq1+1−expAqΔtAqΔt,B1,q=CqexpAqΔtAq1−expAqΔtAqΔt+expAqΔt.From equation (21), it is clear that the SARC approach has a high efficiency in calculating the near-field radiation. This is because this approach avoids the requirement for storing and processing of the complete history of the equivalent source on the output boundary of the FDTD region, and its computational complexity is reduced as only ONmax. Beyond that, the SARC approach has a better accuracy than the general recursive approach, due to the analytical calculation of the convolution integral in the discrete time interval.


## 2.1. Theory of Near-Field Transformation in 2-D Case

Assume thatE and H are respectively the electric and magnetic fields on the output boundary surface of the FDTD domain. According to the equivalence principle, the surface electric and magnetic currents J,M and charges ρs,ρm can be obtained by(1a)J=n^×H,M=−n^×E,(1b)n^⋅μH=ρs,n^⋅μE=ρm,where n^ denotes the unit normal vector on the output boundary. For the 2-D problem, the vector and scalar potentials can be expressed with the 2-D time-domain Green’s function as(2a)Ar,t=μ4π∫Ldl′∫0t−R/cJr′,t′t−t′2−R2/c2dt′,(2b)Fr,t=ε4π∫Ldl′∫0t−R/cMr′,t′t−t′2−R2/c2dt′,(2c)φr,t=14πμ∫Ldl′∫0t−R/cρsr′,t′t−t′2−R2/c2dt′,(2d)φmr,t=14πμ∫Ldl′∫0t−R/cρmr′,t′t−t′2−R2/c2dt′.Here,μ and ε represent the permeability and permittivity, R=r−r′ is the distance from the observation point r to source field r′, and c is the speed of the EM wave. In accordance with the vector and scalar potentials, the 2-D radiation field exterior to the FDTD region can be written as(3a)Er,t=−1ε∇×Fr,t−∇φ−∂∂tAr,t,(3b)Hr,t=1μ∇×Ar,t−∇φm−∂Fr,t∂t.Substituting equation (2) into equation (3) and shifting the time reference backward by R/c, the formals for calculating the 2-D radiation field of the equivalent sources can be obtained from [21], in which, for the TM case, are(4a)Ezr,t+Rc=12π∫Ldl′−μ∫0t∂Jzr′,t′/∂t′t+R/c−t′2−R2/c2dt′−y−y′R2∫0tt+R/c−t′∂Mxr′,t′/∂t′t+R/c−t′2−R2/c2dt′+x−x′R2∫0tt+R/c−t′∂Myr′,t′/∂t′t+R/c−t′2−R2/c2dt′,(4b)Hxr,t+Rc=12π∫Ldl′−ε∫0t∂Mxr′,t′/∂t′t+R/c−t′2−R2/c2dt′−x−x′μR2∫0tt−t′∂ρmr′,t′/∂t′t+R/c−t′2−R2/c2dt′+y−y′R2∫0tt+R/c−t′∂Jzr′,t′/∂t′t+R/c−t′2−R2/c2dt′,(4c)Hyr,t+Rc=12π∫Ldl′−ε∫0t∂Jzr′,t′/∂t′t+R/c−t′2−R2/c2dt′−y−y′μR2∫0tt+R/c−t′∂ρmr′,t′/∂t′t+R/c−t′2−R2/c2dt′+x−x′R2∫0tt+R/c−t′∂Jzr′,t′/∂t′t+R/c−t′2−R2/c2dt′.The corresponding formulas in the TE mode can be derived with the principle of duality, which are not described here. It is obvious that the convolution kernel in equation (4) has a square-root singularity at t′=t. An efficient way in addressing the singularity is to approximately calculate the convolution integral in the interval t−Δt,t using the linear interpolation. Assuming t−t′=τ, the time derivative of the equivalent source, e.g., Mr,t′, in the time-reversed and offset form can be given by(5)MnΔt−τ≈Mn+Mn−1−MnΔtττ∈0,Δt.For instance, substituting (5) into the second term of (4a), we have(6)∫t−Δttt+R/c−t′Mr′,tt+R/c−t′2−R2/c2dt′=Mn+RMn−Mn−1cΔtτ+Rc2−R2c2|τ=0τ=Δt+Mn−1−MnΔtτ+R/c2τ+Rc2−R2c2+R22c2lnτ+Rc+τ+Rc2−R2c2τ=0τ=Δt.A similar procedure can be implemented for other terms in equation (4) to obtain the corresponding approximate results in the interval t−Δt,t. The convolution integral in the interval 0,t−Δt is computed by using the high-efficient SARC approach, as will be elaborated later.


## 2.2. Prony Approximation of Convolution Kernel

Although equation (4) is applicable for the computer execution, the direct integration of it is such cumbersome as to be impractical, due to the need for storing and processing all past values of the time derivative of the equivalent source. The computational complexity of this scheme scales as ONmax2, where Nmax is the total time step. In this paper, the SARC approach will be applied to decrease the complexity of the convolution integral. It can be found from equation (4) that the convolution kernel mainly consists of two forms: 1/t+R/c2−R2/c2 and t+R/c/t+R/c2−R2/c2, which make it difficult to execute the recursive operation of the convolution. Based on the idea of the signal estimation in the high-speed circuit [22], the Prony’s method is used to approximate the convolution kernel with a sum of the exponential terms:(7)fN≈∑i=1QCiexpAiN=∑i=1QCiμiN.Here,fN is the discrete sample of the convolution kernel at t=NΔt, where N denotes the equally spaced time step, Q represents the total number of poles, and Ai and Ci represent the coefficient and amplitude of the ith pole. The Ai and Ci are evaluated in the Prony’s method by solving two sets of linear equations with an intermediate nonlinear equation with Qth degree. Firstly, Nm equally spaced time samples of the convolution kernel are taken to set up a set of linear equation as(8)f0f1f2⋯fQ−1f1f2f3⋯fQ⋯⋯⋯⋱⋯fNm−Q−1fNm−QfNm−Q+1⋯fNm−2α1α2⋯αQ=−fQ−fQ+1⋯−fNm−1,where Nm−Q>Q. The roots α1,α2,…,αQ can be determined by solving this oversampled set of linear equations with the least-squares approach. By virtue of equation (8), the μ can be found as the roots of the following Qth degree equation:(9)μQ+αQμQ−1+αQ−1μQ−2+⋯+α2μ+α1=0.If the complex conjugate roots arise in the solution to equation (9), use their magnitude to substitute for the complex pair of μ, and cancel any duplicate. Substituting μi to equation (7) and introducing Nm discrete samples of the convolution kernel, a new set of the linear equations with the Q values of the Ci can be established as follows:(10)111⋯1μ1μ2μ3⋯μQμ12μ22μ32⋯μQ2⋯⋯⋯⋱⋯μ1Nm−1μ2Nm−1μ3Nm−1⋯μQNm−1C1C2C3⋯CQ=f0f1f2⋯fNm−1,where Nm>Q. The coefficients C1,C2,…,CQ can be obtained by using the least-squares approach to equation (10).Figure2 shows the comparison of the fitting results by the Prony’s approximation with the analytical results, in which the red short dash denotes the analytical result and blue line denotes the fitting result with Prony’s method. For the fitting calculation, the space increment is Δ=0.05m, the time increment is Δt=Δ/2c, and the observation point and source point are located at 20Δ,80Δ and 30Δ,10Δ, respectively. The sampling of the convolution kernel begins at intR/cΔt+1, and the interval is 20Δt. In Figures 2(a) and 2(b), the 1/τ+R/c2−R2/c2 and τ+R/c/τ+R/c2−R2/c2 are respectively expressed as convolution kernels 1 and 2 for convenience. The total numbers of the samples and poles are selected to be 30 and 10, respectively. It can be observed from Figure 2 that the nearly identical results are produced with the two methods, implying the reliability of the Prony’s method.Figure 2
Comparison of the fitting results by the Prony’s approximation with the analytical results. (a) Convolution kernel 1. (b) Convolution kernel 2.
(a)(b)


## 2.3. Semianalytical Recursive Convolution Approach

The RASC approach is a fruitful scheme for the calculation of the convolution integral, which is based on the idea of the efficient analysis of the transfer response in the digital signal processing (DSP) [23]. For a linear time-invariance and dynamical system, as shown in Figure 3, the system response yt can be expressed as the convolution of the input signal xt with the impulse response ht. We have(11)yt=ht∗xt=∫0tht−t′xt′dt′.Figure 3
Transfer response for linear time-invariance system.In general, the impulse response in an asymptotically stable system is of the form,(12)ht=∑q=1QHqexp−βqtUt,where q stands for the number of poles, and Ut is the switch function:(13)Ut=1t≥00t<0.Substituting equation (12) into equation (11), the system response at t=nΔt can be obtained as(14)yn=∑q=1Q∫0nΔtHqexp−βqnΔt−t′xt′dt′.The system response for theqth pole is(15)yqn=∫0nΔtHqexp−βqnΔt−t′xt′dt′.Dividing the integral of equation (15) into two intervals, 0,n−1Δt and n−1Δt,nΔt, we get to the recursive formula of system response as(16)yqn=exp−βqΔtyqn−1+Hqexp−βqnΔt∫n−1ΔtnΔtexpβqt′xt′dt′.In this paper, the near-field radiated by the equivalent source is analogous to the system response in the DSP. The convolution kernel can be regard as the impulse response after it is approximated with a sum of exponential terms. The time derivative of the equivalent source corresponds to the input signal. Thus, the near-field radiation can be updated by applying the SARC approach. Next, focus on the SARC approach for evaluating the radiation field in equation (4), where the representative convolutions are of the forms,(17a)It+Rc=∫0t−ΔtWr′,t′t+R/c−t′2−R2/c2dt′,(17b)Kt+R/c=∫0t−Δtt+R/c−t′Wr′,t′t+R/c−t′2−R2/c2dt′.For equation (17a) as example, using the Prony’s method mentioned above can approximate the convolution kernel as(18)1t+R/c2−R2/c2≈∑q=1QCqexpAqt.According to equation (16), the temporal convolution in the interval 0,n−1Δt can be expressed as(19)Iqn−1=expAqΔtIqn−2+CqexpAqnΔt∫n−2Δtn−1Δtexp−Aqt′Wt′dt′,where Wt can be approximated in accordance with the linear interpolation approach as follows:(20)Wt′≈Wn−1+Wn−1−Wn−2Δtt′−n−2Δt.Substituting equation (20) into equation (19) and analytically solving the integral over the interval n−2Δt,n−1Δt, we find the following field-update equation:(21)Iqn−1=expAqΔtIqn−2+B0,qWn−1+B1,qWn−2,where(22)B0,q=−CqexpAqΔtAq1+1−expAqΔtAqΔt,B1,q=CqexpAqΔtAq1−expAqΔtAqΔt+expAqΔt.From equation (21), it is clear that the SARC approach has a high efficiency in calculating the near-field radiation. This is because this approach avoids the requirement for storing and processing of the complete history of the equivalent source on the output boundary of the FDTD region, and its computational complexity is reduced as only ONmax. Beyond that, the SARC approach has a better accuracy than the general recursive approach, due to the analytical calculation of the convolution integral in the discrete time interval.


## 3. Results and Discussion

Figure4 shows the geometry of the near-field scattering from rough surface. The observation point P is located above the rough surface, and θi denotes the incident angle. Without loss of generality, the one-dimensional Weierestrass-Mandelbrot function is used here to generate the rough surface [24]. In the numerical calculation, the length of the rough surface is taken to be L=160m and the spatial increment is taken to be Δ=Δx=Δy=0.05m to allow for reasonable resolution of the rough surface model discretized with the staircase approximation, and the time increment is set to be Δt=Δ/2c in accordance with Courant stability condition [22]. The 2-D FDTD space lattice is terminated by the uniaxial perfectly matched layer (UPML) [25], and its thickness is chosen to be 10Δ. Consider the rough surface is exposed to the Gaussian pulse, whose expression is Ei=exp−4πt−t0/τ2. The parameters of the incident wave are adopted as follows: t0=60Δt and τ=0.8t0.Figure 4
Geometry of near-field scattering from rough surface.To illustrate the accuracy and efficiency of the proposed hybrid algorithm, Figures5 and 6 present the comparison of the near-field scattering results by the hybrid algorithm with those by the traditional FDTD method for the TE and TM polarized conditions, respectively. In Figures 5(a) and 5(b), the near-zone observation points are respectively located at 0.5m,5m and 0.5m,16m, the incident angle is selected to be θi=200, and the root-mean-square (RMS) height is δ=0.1m, the fractal dimension is D=1.5, and the permittivity of the rough surface is ε=10.801,0.384i. In Figures 6(a) and 6(b), the near-zone observation points are positioned at 1.5m,5m and 1.5m,16m; the rest of the parameters are shown as follows: θi=300, δ=0.2m, D=1.5, and ε=10.801,0.384i. It can be observed from Figures 5 and 6 that the time-domain scattered waveforms for the two methods are in good agreement for different conditions. This result demonstrates the accuracy of the hybrid algorithm proposed in this paper. Beyond that, hybrid algorithm can dramatically lessen the computation cost compared with the traditional FDTD method. For the TE polarized wave incidence and 4000 time-steps updating as example, the computation costs required for two methods are provided in Table 1. When the height of the near-zone observation point is h=5m, the number of FDTD meshes, the memory usage, and the run time in the hybrid algorithm are reduced to about 72%, 76%, and 74% of those in the traditional FDTD method. For h=16m, the key factors above using the hybrid algorithm are only 39%, 42%, and 33% of those using the traditional method, respectively. It is obvious that the hybrid algorithm becomes more efficient with the increase in the height of the near-zone observation point. This is because that the FDTD region in the traditional method has to be greatly enlarged to enclose the rough surface and observation point, whereas that in the hybrid algorithm remains unchanged when the observation point rises. In addition, we can find from Table 1 that the proposed algorithm differs little from the traditional FDTD method in the run time for h=2.5m. This is owing to the fact that, in this case of very low observation point, the advantage of the hybrid algorithm in reducing the FDTD meshes is almost offset by its calculation of the near-field radiation. As the height of the near-zone observation point increases, the proposed hybrid algorithm becomes more favorable in the computation efficiency.Figure 5
Near-fieldHz scattered from the rough surface by FDTD and by hybrid algorithms for the TE case. (a) h=5m; (b) h=16m.
(a)(b)Figure 6
Near-fieldEz scattered from the rough surface by FDTD and by hybrid algorithms for the TM case. (a) h=5m; (b) h=16m.
(a)(b)Table 1
Computation costs of the hybrid algorithm and the FDTD method for the TE case.
FDTDHybrid algorithmFDTDHybrid algorithmFDTDHybrid algorithmObservation point attitudeh=2.5mh=2.5mh=5mh=5mh=16mh=16mNumber of cells6880006080008480006080001552000608000Memory usage (MB)27.925.433.425.460.525.4Run time (s)5165116915121551512∗Calculated by a computer with a 2.30 GHz processor (Intel (R) Xeon (R) Gold 5118 CPU), 64 GB memory.


## 4. Conclusions

In this paper, a novel hybrid algorithm has been proposed to calculate the transient near-field scattering from rough surface for both TE and TM polarization. This hybrid algorithm is based on the combination of the FDTD and the 2-D time-domain Huygens’ principle, and it is the development of the traditional FDTD method. In the hybrid algorithm, the FDTD is only performed in the computational domain containing the rough surface, which makes a large number of FDTD meshes eliminated in the space between the rough surface and observation point. The proposed hybrid algorithm is validated by comparing the near-field scattered waveforms of the rough surface for both the hybrid algorithm and the traditional FDTD method. The numerical results show that the hybrid algorithm is an accurate and efficient technique in the calculation of the transient near-field scattering from the rough surface. Moreover, the relative efficiency of the hybrid algorithm becomes higher with the increase in the height of the observation point. In the future, the hybrid algorithm based on SARC approach will be developed to study the transient composite scattering and coupling of the dispersive target and background.



---

*Source: 7164286-2021-08-26.xml*