Multi-parameter identification and shape reconstruction for unbounded fractal rough surfaces with tapered wave incidence

Abstract

In this paper, we study the inverse scattering of tapered waves from sound soft surfaces or horizontally polarized electromagnetic waves from perfectly conducting surfaces. By the related Fréchet derivative with respect to a given surface of the solution of direct scattering problems, an efficient algorithm is proposed to reconstruct the parameters of unbounded rough surfaces from a set of scattered field measurements. Tapered wave is introduced to realize the asymptotic truncation for unbounded fractal rough surface. In order to improve the results of reconstruction, multi-angle and multi-frequency incidence strategy has been used. The numerical results show that the method yields satisfactory reconstructions for some rough surface profiles.

Keywords:

• Inverse scattering
• unbounded rough surface
• Fréchet derivative
• BFGS algorithm

AMS Subject Classifications:

• 78A46
• 65N21
• 45Q05
• 35J05

1. Introduction

Inverse scattering problems are of great interest because of their potential applications in biomedical imaging, non-destructive testing and diagnostics, remote sensing applications, applied geophysics and non-invasive subsurface monitoring. The inverse scattering of rough surface is such problem that the incident wave is given and the measurement data of scattering wave are assumed to be known, on the basis of which the characteristic information of the unbounded rough surfaces such as geometric parameters and physical properties can be confirmed and reconstructed. Parameters inversion of rough surface constitutes an important class of problems in inverse scattering theory due to the large range of applications such as microwave remote sensing, underwater acoustics and target recognition. In these problems, one tries to recover the multi-parameter of an unknown fractal rough surface from scattered field measurements in a certain domain. It is a simplification of ground (Surface models) imaging via ground-penetrating radar (GPR). In GPR systems, arrays of above-ground transmitters and receivers illuminate surfaces of interest and receive scattered data from the rough surface. The rough surface to be reconstructed can be either perfectly conducting or a dielectric interface separating two dielectric media.

The problem of parameter identification and shape reconstruction has attracted the attention of many scholars. Some investigations have been developed to solve the problems in recent years. For the periodic rough surfaces (such as grating), some scattering and inverse scattering problems have been studied (see, e.g. [1–5] and references therein). In [6], the authors discuss the reconstruction of a fractal rough surface with periodic boundary conditions, a single-parameter (fractal dimension) is reconstructed directly by linear search method. For the locally rough surface, the integral equation method for scattering and inverse scattering problems has been considered (see, e.g. [7,8]). In reality, the ground is an unbounded randomly rough surface, which makes it more challenging to reconstruct its profile. Near-field imaging of one kind of unbounded rough surfaces is considered in [9]. Several analytic and numerical techniques have been developed for studying electromagnetic scattering and inverse scattering from random rough surfaces without and with certain targets, respectively (see, e.g. [10–15]). Many works are based on the Kirchhoff approximation with the rough surface assumed to be locally planar (see, e.g. [16–18]), which accounts only for the single scattering and not the multiple scattering. The other methods such as the coupled integral equation method (see, e.g. [19,20]), the linear sampling method (see, e.g. [21]) and the point source method (see, e.g. [22]) have been proposed for solving the non-linear inverse problem and determining the profile of the rough surface.

In this paper, the rough surface is characterized by one-dimensional band-limited Weierstrass–Mandelbrot function, and it is subject to the Dirichlet (zero field) boundary condition. With fast and accurate solver for the direct scattering problem which we presented in [23], the optimization algorithms for multi-parameters inversion and shape reconstruction of unbounded non-periodic rough surface are investigated. The reconstruction algorithm is based on merging a fast direct solver and a non-linear optimal technique. The goal is to minimize an objective function relating parameters of the rough surface to the scattered field data. Firstly, we introduced the related Fréchet derivative with respect to parameters of a given surface for the solution of direct scattering problems. It is a basis for the optimization algorithm. Secondly, we use BFGS method [24] to solve the non-linear optimization problem. For each step, it needs to apply a direct solver several times for solving the direct problem and the Fréchet derivative of objective function. Finally, the multi-angle and multi-frequency incidence strategy has also been considered. By Abbe’s resolution formula, it is reasonable to expect the better results should be obtained when the multi-frequency data are used. Additionally, illuminating the unknown rough surface from more directions should make the surface easier to recover. The measurements of scattered field are generated by certain numerical method (such as the regularized conjugate gradient method (RCGM),[23] and method of moments [15]), and they are used to replace the real collected GPR-type data to test the reconstruction algorithm. Several key issues will be examined such as the reconstruction accuracy, the influence of measurement error, initial data and the surface roughness.

The organization of the paper is as follows: In Section 2, we describe the scattering problem of unbounded rough surface with tapered wave incidence and the boundary integral formalism. To solve the optimization problem for the parameters inversion of the function f from these measured data, the Fréchet derivative of the scattered field is introduced. The numerical solution of the integral equation is also briefly discussed. In Section 3, the numerical inversion method for reconstructing the rough surface is presented. The feasibility of the method is illustrated through some numerical examples, and the numerical results are discussed in Section 4. Finally, the concluding remarks are given in Section 5.

2. Scattering problem of fractal rough surfaces

2.1. Rough surface function and Scattering problem

The parameters inversion and shape reconstruction of the unbounded rough surface begin with assuming a mathematical model of the surface that can approximately define its height variation. The geometric shape of unbounded rough surfaces (such as sea surface and ground surface) are often represented by means of a deterministic function as well as a stochastic process (see, e.g. [14,25]). Since fractals hold in balance long-range order and short-range disorder, here the band-limited fractal function f is used to describe the rough surface(1) $$\begin{array}{c}\hfill f(x)=\widehat{h}C\sum _{n=0}^{M_f-1}{b}^{(\widehat{D}-2)n}sin(K{b}^{n}x+{\mathit{\phi }}_n),x\in \mathbb{R},\end{array}$$(1)

where $\mathrm{\Lambda }$ is the fundamental surface period, $K=2\mathit{\pi }/\mathrm{\Lambda }$ is the fundamental surface wave number, $b(b>1)$ is the frequency scaling parameter, $\widehat{h}$ is the root mean square height, $\widehat{D}(1<\widehat{D}<2)$ is the fractal dimension, $M_f$ is the highest harmonic wave number, $C=\sqrt{2(1-b^{2(\widehat{D}-2)})/(1-b^{2M_f(\widehat{D}-2)})}$ is a normalization factor and ${\mathit{\phi }}_n$ is the specified phases. For more details of the band-limited fractal function, we refer the interested reader to [6,26,27].

We assume that the unbounded fractal rough surfaces $\mathrm{\Gamma }$ are impenetrable and $\mathrm{\Omega }$ is the propagation domain above f, which can be expressed as$$\mathrm{\Gamma }=\{(x,z)\in {\mathbb{R}}^2|z=f(x)\},\mathrm{\Omega }=\{(x,z)\in {\mathbb{R}}^2|z>f(x)\}.$$

For considering the scattering problem of unbounded rough surfaces, the Thorsos tapered wave ${\mathrm{\Psi }}^{inc}(x,z)$ is used as the incident field(2) $$\begin{array}{c}\hfill {\mathrm{\Psi }}^{inc}(x,z)=exp(\mathrm{i}k(xsin{\mathit{\theta }}_{inc}-zcos{\mathit{\theta }}_{inc})(1+w(x,z)))exp(\frac{-{(x+ztan{\mathit{\theta }}_{inc})}^2}{g^2}),\end{array}$$(2)

which satisfies the non-homogeneous Helmholtz equation(3) $$\begin{array}{c}\hfill \mathrm{\Delta }{\mathrm{\Psi }}^{inc}+k^2{\mathrm{\Psi }}^{inc}=k^{2}R\mathrm{in}\mathrm{\Omega },\end{array}$$(3)

where(4) $$\begin{array}{ccc}\hfill R\equiv {\mathrm{\Psi }}^{inc}(x,z)\{& -\hfill & \hfill w^2(\overline{r})-16\frac{{(xsin{\mathit{\theta }}_{inc}-zcos{\mathit{\theta }}_{inc})}^2{(x+ztan{\mathit{\theta }}_{inc})}^2}{k^4{g}^8{cos}^6{\mathit{\theta }}_{inc}}\\ \hfill & +\frac{4\mathrm{i}k(xsin{\mathit{\theta }}_{inc}-zcos{\mathit{\theta }}_{inc})}{k^4{g}^4{cos}^4{\mathit{\theta }}_{inc}}(1-\frac{4{(x+ztan{\mathit{\theta }}_{inc})}^2}{g^2})\}\hfill \end{array}$$(4)

and(5) $$\begin{array}{c}\hfill w(x,z)={(kgcos{\mathit{\theta }}_{inc})}^{-2}(2{(x+ztan{\mathit{\theta }}_{inc})}^2{g}^{-2}-1).\end{array}$$(5)

In fact, $|R/{\mathrm{\Psi }}^{inc}|=O(1/k^3{g}^3{cos}^3{\mathit{\theta }}_{inc})$ as $kg\to \infty$. ${\mathit{\theta }}_{inc}(0<|{\mathit{\theta }}_{inc}|<\frac{\mathit{\pi }}{2})$ is the angle of incidence (it is the angle between the direction of propagation and the negative z axis), and g is the parameter that controls the tapering, $k=2\mathit{\pi }/\mathit{\lambda }$ is the wave number, $\mathit{\lambda }$ is the wavelength of incident wave, $kgcos{\mathit{\theta }}_{inc}\gg 1$ and $\mathrm{i}=\sqrt{-1}$. The tapered wave is introduced by Thorsos in 1988 [16], which has been well received in the electromagnetics community. In the process of numerical calculation, the rough surface must be truncated in a finite interval. This means that the surface current is forced to be zero outside the interval for plane wave incident. If there is an abrupt change of surface current from non-zero to zero, artificial reflection from the two endpoints will occur. Note that at $z=0$ (or z is equal to a constant), the amplitude of the tapered incident field decreased exponentially for $|x|>g$. So it can be used to realize the asymptotic truncation.

The total field $\mathrm{\Psi }$ consists of the incident field ${\mathrm{\Psi }}^{inc}$ and the scattered field ${\mathrm{\Psi }}^s$:(6) $$\begin{array}{c}\hfill \mathrm{\Psi }={\mathrm{\Psi }}^{inc}+{\mathrm{\Psi }}^s,\end{array}$$(6)

where the scattered field ${\mathrm{\Psi }}^s$ satisfies the homogeneous Helmholtz equation(7) $$\begin{array}{c}\hfill \mathrm{\Delta }{\mathrm{\Psi }}^s+k^2{\mathrm{\Psi }}^s=0\mathrm{in}\mathrm{\Omega }.\end{array}$$(7)

We describe the scattering model mathematically as follows (see, e.g. [16,23,28,29])(8) $$\begin{array}{c}\hfill \mathrm{\Delta }\mathrm{\Psi }+k^2\mathrm{\Psi }=k^{2}R\mathrm{in}\mathrm{\Omega }.\end{array}$$(8)

For a sound soft (or perfectly conducting) surface, the total filed $\mathrm{\Psi }$ vanishes on the rough surface,(9) $$\begin{array}{c}\hfill {\mathrm{\Psi }|}_{\mathrm{\Gamma }}=0.\end{array}$$(9)

The scattering problem is not complete without a suitable radiation condition. Physically, such a condition ensures that the scattered field is propagating away from the surfaces while mathematically, it ensures uniqueness of the solution. A natural approach to obtain a radiation condition, often used in applications, is to consider the field ${\mathrm{\Psi }}^s$ only in the domain ${\mathrm{\Omega }}_l:=\{(x,z)\in {\mathbb{R}}^2|z\ge l>max\{0,f(x)\}\}$ and to compute its Fourier transform $(\mathcal{F}{\mathrm{\Psi }}^s)(\mathit{\xi },z):=\frac{1}{\sqrt{2\mathit{\pi }}}{\int }_{-\infty }^{+\infty }exp(-\mathrm{i}\mathit{\xi }x){\mathrm{\Psi }}^s(\mathit{\xi },z)\mathrm{d}\mathit{\xi }$ with respect to the x variable. This leads to an ordinary differential equation(10) $$\begin{array}{c}\hfill \frac{\mathit{\partial }\mathcal{F}{\mathrm{\Psi }}^s}{\mathit{\partial }{z}^2}(\mathit{\xi },z)+(k^2-{\mathit{\xi }}^2)\mathcal{F}{\mathrm{\Psi }}^s(\mathit{\xi },z)=0,\mathit{\xi }\in \mathbb{R},\end{array}$$(10)

with some initial condition, e.g. $\mathcal{F}{\mathrm{\Psi }}_{Dir}(\mathit{\xi })=(\mathcal{F}{\mathrm{\Psi }}^s)(\mathit{\xi },l)\in L^2(\mathbb{R})$ for ${\mathrm{\Psi }}_{Dir}(\mathit{\xi })={\mathrm{\Psi }}^s(\mathit{\xi },l),x\in \mathbb{R}$. Selecting the outgoing of the two linearly independent solutions, we arrive at(11) $$\begin{array}{c}\hfill {\mathrm{\Psi }}^s(x,z)=\frac{1}{\sqrt{2\mathit{\pi }}}{\int }_{-\infty }^{+\infty }exp(\mathrm{i}x\mathit{\xi }+\mathrm{i}(z-l)\sqrt{k^2-{\mathit{\xi }}^2})(\mathcal{F}{\mathrm{\Psi }}_{Dir})(\mathit{\xi })\mathrm{d}\mathit{\xi }in{\mathrm{\Omega }}_l.\end{array}$$(11)

Equation (11(11) $$\begin{array}{c}\hfill {\mathrm{\Psi }}^s(x,z)=\frac{1}{\sqrt{2\mathit{\pi }}}{\int }_{-\infty }^{+\infty }exp(\mathrm{i}x\mathit{\xi }+\mathrm{i}(z-l)\sqrt{k^2-{\mathit{\xi }}^2})(\mathcal{F}{\mathrm{\Psi }}_{Dir})(\mathit{\xi })\mathrm{d}\mathit{\xi }in{\mathrm{\Omega }}_l.\end{array}$$(11) ) is often referred to as an angular spectrum representation. The formalism can easily be made mathematically rigorous if $\mathcal{F}{\mathrm{\Psi }}_{Dir}(\mathit{\xi })\in L^2(\mathbb{R})$. It is worth noting that the radiation conditions (11(11) $$\begin{array}{c}\hfill {\mathrm{\Psi }}^s(x,z)=\frac{1}{\sqrt{2\mathit{\pi }}}{\int }_{-\infty }^{+\infty }exp(\mathrm{i}x\mathit{\xi }+\mathrm{i}(z-l)\sqrt{k^2-{\mathit{\xi }}^2})(\mathcal{F}{\mathrm{\Psi }}_{Dir})(\mathit{\xi })\mathrm{d}\mathit{\xi }in{\mathrm{\Omega }}_l.\end{array}$$(11) ) is different from the classical Sommerfeld’s radiation conditions. The Sommerfeld’s radiation condition is not a valid characterization of outgoing waves for scattering problems at unbounded rough surfaces. A simple example is the scattering on a flat surface: for a downward propagating incident wave, the scattered field ${\mathrm{\Psi }}^s(x,z)$ is simply some reflected wave and hence satisfies Sommerfeld’s radiation condition only in the propagation direction. So the angular spectrum representation is proposed to replace the Sommerfeld’s radiation condition (see, e.g. [30–32]).

2.2. Boundary integral equations

The boundary integral equation (BIE) approach has been developed to handle the scattering problem (8(8) $$\begin{array}{c}\hfill \mathrm{\Delta }\mathrm{\Psi }+k^2\mathrm{\Psi }=k^{2}R\mathrm{in}\mathrm{\Omega }.\end{array}$$(8) ), (9(9) $$\begin{array}{c}\hfill {\mathrm{\Psi }|}_{\mathrm{\Gamma }}=0.\end{array}$$(9) ) and (11(11) $$\begin{array}{c}\hfill {\mathrm{\Psi }}^s(x,z)=\frac{1}{\sqrt{2\mathit{\pi }}}{\int }_{-\infty }^{+\infty }exp(\mathrm{i}x\mathit{\xi }+\mathrm{i}(z-l)\sqrt{k^2-{\mathit{\xi }}^2})(\mathcal{F}{\mathrm{\Psi }}_{Dir})(\mathit{\xi })\mathrm{d}\mathit{\xi }in{\mathrm{\Omega }}_l.\end{array}$$(11) ). Using Green’s theorem and the angular spectrum representation radiation condition, we derived the BIE (12(12) $$\begin{array}{cc}& {\mathrm{\Psi }}^{inc}(\overline{r})+{\int }_{\mathrm{\Gamma }}(\frac{\mathit{\partial }\mathrm{\Phi }(\overline{r},\overline{r^{\prime }})}{\mathit{\partial }{n}^{\prime }}\mathrm{\Psi }(\overline{r^{\prime }})-\mathrm{\Phi }(\overline{r},\overline{r^{\prime }})\frac{\mathit{\partial }\mathrm{\Psi }(\overline{r^{\prime }})}{\mathit{\partial }{n}^{\prime }})\mathrm{d}{s}^{\prime }+{\int }_{\mathrm{\Omega }}{k}^{2}R(\overline{r^{\prime }})\mathrm{\Phi }(\overline{r},\overline{r^{\prime }})\mathrm{d}\mathit{\sigma }\hfill \\ \hfill & =\left\{\begin{array}{cc}\mathrm{\Psi }(\overline{r}),\hfill & \overline{r}\in \mathrm{\Omega },\hfill \\ 0,\hfill & \overline{r}\in {\mathbb{R}}^2\backslash \mathrm{\Omega },\hfill \end{array}\right.\hfill \end{array}$$(12) ) in [33](12) $$\begin{array}{cc}& {\mathrm{\Psi }}^{inc}(\overline{r})+{\int }_{\mathrm{\Gamma }}(\frac{\mathit{\partial }\mathrm{\Phi }(\overline{r},\overline{r^{\prime }})}{\mathit{\partial }{n}^{\prime }}\mathrm{\Psi }(\overline{r^{\prime }})-\mathrm{\Phi }(\overline{r},\overline{r^{\prime }})\frac{\mathit{\partial }\mathrm{\Psi }(\overline{r^{\prime }})}{\mathit{\partial }{n}^{\prime }})\mathrm{d}{s}^{\prime }+{\int }_{\mathrm{\Omega }}{k}^{2}R(\overline{r^{\prime }})\mathrm{\Phi }(\overline{r},\overline{r^{\prime }})\mathrm{d}\mathit{\sigma }\hfill \\ \hfill & =\left\{\begin{array}{cc}\mathrm{\Psi }(\overline{r}),\hfill & \overline{r}\in \mathrm{\Omega },\hfill \\ 0,\hfill & \overline{r}\in {\mathbb{R}}^2\backslash \mathrm{\Omega },\hfill \end{array}\right.\hfill \end{array}$$(12)

where $\mathrm{\Phi }(\overline{r},\overline{r^{\prime }}):=\frac{i}{4}{H}_0^{(1)}(k|\overline{r}-\overline{r^{\prime }}|),$ it is the fundamental solution to the Helmholtz equation in ${\mathbb{R}}^2$, $n^{\prime }=\frac{-f^{\prime }(x^{\prime })\widehat{x}+\widehat{z}}{\sqrt{1+{[f^{\prime }(x^{\prime })]}^2}}$ is the unit normal vector pointing towards the inside of $\mathrm{\Omega }$ from $\mathrm{\Gamma }$, where $\widehat{x}$ and $\widehat{z}$ represent the unit vectors of x axis and z axis, respectively.

It is notable that most of the literature considers the approximation for the tapered wave rough surface scattering problem as follows:(13) $$\begin{array}{c}\hfill {\mathrm{\Psi }}^{inc}(\overline{r})+{\int }_{\mathrm{\Gamma }}(\frac{\mathit{\partial }\mathrm{\Phi }(\overline{r},\overline{r^{\prime }})}{\mathit{\partial }{n}^{\prime }}\mathrm{\Psi }(\overline{r^{\prime }})-\mathrm{\Phi }(\overline{r},\overline{r^{\prime }})\frac{\mathit{\partial }\mathrm{\Psi }(\overline{r^{\prime }})}{\mathit{\partial }{n}^{\prime }})\mathrm{d}{s}^{\prime }=\left\{\begin{array}{cc}\mathrm{\Psi }(\overline{r}),\hfill & \overline{r}\in \mathrm{\Omega },\hfill \\ 0,\hfill & \overline{r}\in {\mathbb{R}}^2\backslash \mathrm{\Omega }.\hfill \end{array}\right.\end{array}$$(13)

The derivation of BIEs for scattering by unbounded rough surface with plane wave incidence has been studied in detail by J.A. DeSanto and P.A. Martin in [31,32,34]. They pointed out that the BIE (13(13) $$\begin{array}{c}\hfill {\mathrm{\Psi }}^{inc}(\overline{r})+{\int }_{\mathrm{\Gamma }}(\frac{\mathit{\partial }\mathrm{\Phi }(\overline{r},\overline{r^{\prime }})}{\mathit{\partial }{n}^{\prime }}\mathrm{\Psi }(\overline{r^{\prime }})-\mathrm{\Phi }(\overline{r},\overline{r^{\prime }})\frac{\mathit{\partial }\mathrm{\Psi }(\overline{r^{\prime }})}{\mathit{\partial }{n}^{\prime }})\mathrm{d}{s}^{\prime }=\left\{\begin{array}{cc}\mathrm{\Psi }(\overline{r}),\hfill & \overline{r}\in \mathrm{\Omega },\hfill \\ 0,\hfill & \overline{r}\in {\mathbb{R}}^2\backslash \mathrm{\Omega }.\hfill \end{array}\right.\end{array}$$(13) ) is an approximate equation with tapered wave incidence (see [31], Sec.VI.D). We consider the integral equation (12(12) $$\begin{array}{cc}& {\mathrm{\Psi }}^{inc}(\overline{r})+{\int }_{\mathrm{\Gamma }}(\frac{\mathit{\partial }\mathrm{\Phi }(\overline{r},\overline{r^{\prime }})}{\mathit{\partial }{n}^{\prime }}\mathrm{\Psi }(\overline{r^{\prime }})-\mathrm{\Phi }(\overline{r},\overline{r^{\prime }})\frac{\mathit{\partial }\mathrm{\Psi }(\overline{r^{\prime }})}{\mathit{\partial }{n}^{\prime }})\mathrm{d}{s}^{\prime }+{\int }_{\mathrm{\Omega }}{k}^{2}R(\overline{r^{\prime }})\mathrm{\Phi }(\overline{r},\overline{r^{\prime }})\mathrm{d}\mathit{\sigma }\hfill \\ \hfill & =\left\{\begin{array}{cc}\mathrm{\Psi }(\overline{r}),\hfill & \overline{r}\in \mathrm{\Omega },\hfill \\ 0,\hfill & \overline{r}\in {\mathbb{R}}^2\backslash \mathrm{\Omega },\hfill \end{array}\right.\hfill \end{array}$$(12) ) directly. With Dirichlet boundary condition, the surface integral Equation (12(12) $$\begin{array}{cc}& {\mathrm{\Psi }}^{inc}(\overline{r})+{\int }_{\mathrm{\Gamma }}(\frac{\mathit{\partial }\mathrm{\Phi }(\overline{r},\overline{r^{\prime }})}{\mathit{\partial }{n}^{\prime }}\mathrm{\Psi }(\overline{r^{\prime }})-\mathrm{\Phi }(\overline{r},\overline{r^{\prime }})\frac{\mathit{\partial }\mathrm{\Psi }(\overline{r^{\prime }})}{\mathit{\partial }{n}^{\prime }})\mathrm{d}{s}^{\prime }+{\int }_{\mathrm{\Omega }}{k}^{2}R(\overline{r^{\prime }})\mathrm{\Phi }(\overline{r},\overline{r^{\prime }})\mathrm{d}\mathit{\sigma }\hfill \\ \hfill & =\left\{\begin{array}{cc}\mathrm{\Psi }(\overline{r}),\hfill & \overline{r}\in \mathrm{\Omega },\hfill \\ 0,\hfill & \overline{r}\in {\mathbb{R}}^2\backslash \mathrm{\Omega },\hfill \end{array}\right.\hfill \end{array}$$(12) ) is converted to(14) $$\begin{array}{c}\hfill {\int }_{\mathrm{\Gamma }}\mathrm{\Phi }(\overline{r},\overline{r^{\prime }})\frac{\mathit{\partial }\mathrm{\Psi }(\overline{r^{\prime }})}{\mathit{\partial }{n}^{\prime }}\mathrm{d}{s}^{\prime }={\mathrm{\Psi }}^{inc}(\overline{r})+{\mathit{\rho }}_{\mathrm{\Gamma }}(\overline{r}),\overline{r}\in \mathrm{\Gamma },\end{array}$$(14)

where ${\mathit{\rho }}_{\mathrm{\Gamma }}(\overline{r})={\int }_{\mathrm{\Omega }}{k}^{2}R(\overline{r^{\prime }})\mathrm{\Phi }(\overline{r},\overline{r^{\prime }})\mathrm{d}\mathit{\sigma }$. It remains to discuss the discretization and solution of (14(14) $$\begin{array}{c}\hfill {\int }_{\mathrm{\Gamma }}\mathrm{\Phi }(\overline{r},\overline{r^{\prime }})\frac{\mathit{\partial }\mathrm{\Psi }(\overline{r^{\prime }})}{\mathit{\partial }{n}^{\prime }}\mathrm{d}{s}^{\prime }={\mathrm{\Psi }}^{inc}(\overline{r})+{\mathit{\rho }}_{\mathrm{\Gamma }}(\overline{r}),\overline{r}\in \mathrm{\Gamma },\end{array}$$(14) ). For this, we assume the boundary $\mathrm{\Gamma }$ is parameterized, and it is truncated in an finite integration interval $[-L,L]$(15) $$\begin{array}{c}\hfill {\int }_{-L}^{+L}\mathrm{\Phi }(x,f(x);x^{\prime },f(x^{\prime }))\sqrt{1+{[f^{\prime }(x^{\prime })]}^2}[\frac{\mathit{\partial }\mathrm{\Psi }(\overline{r^{\prime }})}{\mathit{\partial }{n}^{\prime }}{]}_{z^{\prime }=f(x^{\prime })}\mathrm{d}{x}^{\prime }=g^{\mathit{\delta }}(x),x\in [-L,L],\end{array}$$(15)

where $g^{\mathit{\delta }}(x)$ is the corresponding right-hand side. We discretize the integral equation by a collocation method, with N equispaced points corresponding to $x_j=-L+2(j-1)L/N$ on the boundary for $j=1,\dots ,N$. We assign an unknown value $\frac{\mathit{\partial }\mathrm{\Psi }}{\mathit{\partial }n}(x_j)$ at each such point and enforce the integral equation pointwise at the same points. Because the integral is logarithmically singular, the small argument approximation of the Hankel function is employed (see, e.g. [15]). The discretization would yield a dense matrix($N\times N$). The discrete linear systems can be expressed as(16) $$\begin{array}{c}\hfill AU=g^{\overline{\mathit{\delta }}},\end{array}$$(16)

where $\Vert g^{\overline{\mathit{\delta }}}{-g\Vert }_{\infty }={max}_{1\le m\le N}|{\mathit{\rho }}_m|$. For solving the equations, one of the limitations is that the computational complexity grows dramatically as the size of the scattering surface increases. The RCGM with fast multipole acceleration which we proposed in [23] is used to solve the problem.

3. Regularized conjugate gradient method

(Algorithm 3.1 in [23])

(1)	Input the noise level $\overline{\mathit{\delta }}$ and the largest admissible number of iteration steps $m_{max}$, chose constant $\mathit{\tau }$$(\mathit{\tau }>1)$, set $m=1$;
(2)	Input the initial guess $U_0^{\overline{\mathit{\delta }}}$, compute $r_0=g^{\overline{\mathit{\delta }}}-A{U}_0^{\overline{\mathit{\delta }}},p_1=s_0=A^{\ast }{r}_0$;
(3)	Compute the received power by the rough surface$$P^{inc}=\frac{cos{\mathit{\theta }}_i}{2{\mathit{\eta }}_0}g\sqrt{\frac{\mathit{\pi }}{2}}(1-\frac{1+2{tan}^2{\mathit{\theta }}_{inc}}{2{(kgcos{\mathit{\theta }}_{inc})}^2});$$
(4)	Input $E_0^{\overline{\mathit{\delta }}}=0$ and the energy stopping tolerance $\mathit{\epsilon }$;
(5)	do while ($E_{m-1}^{\overline{\mathit{\delta }}}<1+\mathit{\epsilon }$)
(6)	if $(m\ge m_{max}$ or $\Vert s_m{\Vert }_{\infty }\le \mathit{\tau }\overline{\mathit{\delta }}<\Vert s_{m-1}{\Vert }_{\infty }),$ break; else compute$$\begin{array}{cc}& q_m=A{p}_m,\hfill \\ \hfill & {\mathit{\alpha }}_m=\Vert s_{m-1}{\Vert }^2/{\Vert q_m\Vert }^2,\hfill \\ \hfill & U_m^{\overline{\mathit{\delta }}}=U_{m-1}^{\overline{\mathit{\delta }}}+{\mathit{\alpha }}_k{p}_m,\hfill \\ \hfill & d_m=d_{m-1}-{\mathit{\alpha }}_m{q}_m,\hfill \\ \hfill & s_m=A^{\ast }{d}_m,\hfill \\ \hfill & {\mathit{\beta }}_m=\Vert s_m{\Vert }^2/{\Vert s_{m-1}\Vert }^2,\hfill \\ \hfill & p_{m+1}=s_m+{\mathit{\beta }}_m{p}_m;\hfill \end{array}$$
(7)	compute$$\begin{array}{cc}\hfill {\mathrm{\Psi }}_{s,m}^{N,\overline{\mathit{\delta }}}({\mathit{\theta }}_s)& ={\int }_{-L}^L-U_m^{\overline{\mathit{\delta }}}(x^{{}^{\prime }})exp(-\mathrm{i}k[sin{\mathit{\theta }}_s{x}^{{}^{\prime }}+cos{\mathit{\theta }}_{s}f(x^{{}^{\prime }})])\mathrm{d}{x}^{{}^{\prime }},\hfill \\ \hfill {\mathit{\sigma }}_m^{\overline{\mathit{\delta }}}({\mathit{\theta }}_s)& =\frac{|{\mathrm{\Psi }}_{s,m}^{N,\overline{\mathit{\delta }}}({\mathit{\theta }}_s){|}^2}{16\mathit{\pi }k{\mathit{\eta }}_0{P}^{inc}},\hfill \\ \hfill E_m^{\overline{\mathit{\delta }}}& ={\int }_{-\frac{\mathit{\pi }}{2}}^{\frac{\mathit{\pi }}{2}}\mathrm{d}{\mathit{\theta }}_s{\mathit{\sigma }}_m^{\overline{\mathit{\delta }}}({\mathit{\theta }}_s);\hfill \end{array}$$
(8)	Set $m=m+1$; end if end do

It should be noted that steps 3, 4, 5 and 7 in the algorithm can ensure the energy relation for scattering calculation. For more details of the algorithm, we refer to [23].

3.1. The Fréchet derivative of the scattered field

Let $\overline{v}={(v_1,v_2,\dots ,v_n)}^T$ denote the unknown parameter vector of the rough surface. ${\mathit{\theta }}_i$ denotes a sufficiently small admissible perturbation of $v_i$ which is a component of the vector $\overline{v}$, and $v_{i,\mathit{\theta }}:=v_i+{\mathit{\theta }}_i(i=1,2,\dots ,n)$. The corresponding solution of scattering problem is denoted as ${\mathrm{\Psi }}_{i,\mathit{\theta }}^s={\mathrm{\Psi }}_{i,\mathit{\theta }}^s(x,z;v_{i,\mathit{\theta }})$ in ${\mathrm{\Omega }}_{i,\mathit{\theta }}$, where ${\mathrm{\Omega }}_{i,\mathit{\theta }}$ is a domain which is determined by $f(x;v_{i,\mathit{\theta }})$, this is$${\mathrm{\Omega }}_{i,\mathit{\theta }}:=\{(x,z)\in {\mathbb{R}}^2|z>f_{i,\mathit{\theta }}(x)=f(x;v_{i,\mathit{\theta }})\},{\mathrm{\Gamma }}_{i,\mathit{\theta }}:=\{(x,z)\in {\mathbb{R}}^2|z=f_{i,\mathit{\theta }}(x)\}.$$

According to the above definition, we have$${\mathrm{\Psi }}_{i,\mathit{\theta }}^s{|}_{{\mathit{\theta }}_i=0}={\mathrm{\Psi }}^s,{\mathrm{\Omega }}_{i,\mathit{\theta }}{{|}_{{\mathit{\theta }}_i=0}=\mathrm{\Omega },{\mathrm{\Gamma }}_{i,\mathit{\theta }}|}_{{\mathit{\theta }}_i=0}=\mathrm{\Gamma },i=1,2,\dots ,n.$$

The direct scattering problem of rough surface can be formulated as the boundary value problem:(17) $$\begin{array}{c}\hfill \left\{\begin{array}{c}\mathrm{\Delta }{\mathrm{\Psi }}_{i,\mathit{\theta }}^s+k^2{\mathrm{\Psi }}_{i,\mathit{\theta }}^s=0\mathrm{in}{\mathrm{\Omega }}_{i,\mathit{\theta }},\hfill \\ {\mathrm{\Psi }}_{i,\mathit{\theta }}^s=-{\mathrm{\Psi }}^{inc}\mathrm{on}{\mathrm{\Gamma }}_{i,\mathit{\theta }},\hfill \\ {\mathrm{\Psi }}_{i,\mathit{\theta }}^s(x,z;v_{i,\mathit{\theta }})=\frac{1}{\sqrt{2\mathit{\pi }}}{\int }_{-\infty }^{+\infty }exp(\mathrm{i}x\mathit{\xi }+\mathrm{i}(z-l_{i,\mathit{\theta }})\sqrt{k^2-{\mathit{\xi }}^2})(\mathcal{F}{\mathrm{\Psi }}_{i,\mathit{\theta }Dir}^s)(\mathit{\xi };v_{i,\mathit{\theta }})\mathrm{d}\mathit{\xi }\mathrm{in}{\mathrm{\Omega }}_{l_{i,\mathit{\theta }}},\hfill \end{array}\right.\end{array}$$(17)

where $\mathcal{F}{\mathrm{\Psi }}_{i,\mathit{\theta }Dir}^s(\mathit{\xi };v_{i,\mathit{\theta }}):=(\mathcal{F}{\mathrm{\Psi }}_{i,\mathit{\theta }}^s)(\mathit{\xi },l_{i,\mathit{\theta }};v_{i,\mathit{\theta }})\in L^2(\mathbb{R})$ is the Fourier transform of ${\mathrm{\Psi }}_{i,\mathit{\theta }Dir}^s(x;v_{i,\mathit{\theta }}):={\mathrm{\Psi }}_{i,\mathit{\theta }}^s(x,l_{i,\mathit{\theta }};v_{i,\mathit{\theta }})$, $x\in \mathbb{R}$, $l_{i,\mathit{\theta }}>max\{0,f_{i,\mathit{\theta }}(x)\}$, ${\mathrm{\Omega }}_{l_{i,\mathit{\theta }}}:=\{\overline{r}=(x,z)\in {\mathbb{R}}^2|z\ge l_{i,\mathit{\theta }}\}$.

We formulate the problem as the operator forms. For a fixed incident field ${\mathrm{\Psi }}^{inc}$, define a scattering operator $\mathcal{K}:{\mathbb{R}}^n\to H^{\mathit{\sigma }}({\mathrm{\Gamma }}_H)$ ($\mathit{\sigma }\ge \frac{1}{2}$) by(18) $$\begin{array}{c}\hfill \mathcal{K}(\overline{v})={\mathrm{\Psi }}^s(x,H),\end{array}$$(18)

where ${\mathrm{\Gamma }}_H:=\{(x,H)\in {\mathbb{R}}^2|H>{max}_{x\in \mathbb{R}}\{0,f(x;\overline{v})\}\}$, for $x\in \mathbb{R}$, ${\mathrm{\Psi }}^s(x,H)$ is the values of the scattered field ${\mathrm{\Psi }}^s$ on the line $z=H$. In the following, we assume that the scattered field is continuously Fréchet differentiable with respect to parameter $v_i(i=1,2,\dots ,n)$. The Fréchet derivative is denoted as ${\mathrm{\Psi }}_i^{\prime }:=\frac{\mathit{\partial }{\mathrm{\Psi }}_{i,\mathit{\theta }}^s}{\mathit{\partial }{v}_{i,\mathit{\theta }}}{|}_{{\mathit{\theta }}_i=0}=\frac{\mathit{\partial }{\mathrm{\Psi }}_{\overline{v}}^s}{\mathit{\partial }{v}_i}$. Then, we have

Lemma 2.1:

The Fréchet derivative ${\mathrm{\Psi }}_i^{\prime }$ satisfies the following equation:$$\mathrm{\Delta }{\mathrm{\Psi }}_i^{\prime }+k^2{\mathrm{\Psi }}_i^{\prime }=0in\mathrm{\Omega }.$$

Proof:

Let(19) $$\begin{array}{c}\hfill A=\mathrm{\Delta }+k^{2}I,\end{array}$$(19)

where I is the identity mapping. And define(20) $$\begin{array}{c}\hfill {\mathit{\varphi }}_{i,\mathit{\theta }}^s(x,z;v_{i,\mathit{\theta }}):=A{\mathrm{\Psi }}_{i,\mathit{\theta }}^s(x,z;v_{i,\mathit{\theta }}),\mathrm{in}{\mathrm{\Omega }}_{i,\mathit{\theta }}.\end{array}$$(20)

It follows from (17(17) $$\begin{array}{c}\hfill \left\{\begin{array}{c}\mathrm{\Delta }{\mathrm{\Psi }}_{i,\mathit{\theta }}^s+k^2{\mathrm{\Psi }}_{i,\mathit{\theta }}^s=0\mathrm{in}{\mathrm{\Omega }}_{i,\mathit{\theta }},\hfill \\ {\mathrm{\Psi }}_{i,\mathit{\theta }}^s=-{\mathrm{\Psi }}^{inc}\mathrm{on}{\mathrm{\Gamma }}_{i,\mathit{\theta }},\hfill \\ {\mathrm{\Psi }}_{i,\mathit{\theta }}^s(x,z;v_{i,\mathit{\theta }})=\frac{1}{\sqrt{2\mathit{\pi }}}{\int }_{-\infty }^{+\infty }exp(\mathrm{i}x\mathit{\xi }+\mathrm{i}(z-l_{i,\mathit{\theta }})\sqrt{k^2-{\mathit{\xi }}^2})(\mathcal{F}{\mathrm{\Psi }}_{i,\mathit{\theta }Dir}^s)(\mathit{\xi };v_{i,\mathit{\theta }})\mathrm{d}\mathit{\xi }\mathrm{in}{\mathrm{\Omega }}_{l_{i,\mathit{\theta }}},\hfill \end{array}\right.\end{array}$$(17) ) that(21) $$\begin{array}{c}\hfill {\mathit{\varphi }}_{i,\mathit{\theta }}^s(x,z;v_{i,\mathit{\theta }})=0,(x,z)\in {\mathrm{\Omega }}_{i,\mathit{\theta }},\end{array}$$(21)

solving the derivative of (21(21) $$\begin{array}{c}\hfill {\mathit{\varphi }}_{i,\mathit{\theta }}^s(x,z;v_{i,\mathit{\theta }})=0,(x,z)\in {\mathrm{\Omega }}_{i,\mathit{\theta }},\end{array}$$(21) ) about ${\mathit{\theta }}_i$ and taking derivative values at ${\mathit{\theta }}_i=0$, we have(22) $$\begin{array}{c}\hfill \frac{\mathit{\partial }{\mathit{\varphi }}_{i,\mathit{\theta }}^s}{\mathit{\partial }{\mathit{\theta }}_i}(x,z;v_{i,\mathit{\theta }}){|}_{{\mathit{\theta }}_i=0}=\frac{\mathit{\partial }{\mathit{\varphi }}_{i,\mathit{\theta }}^s}{\mathit{\partial }z}(x,z;v_{i,\mathit{\theta }})\frac{\mathit{\partial }z}{\mathit{\partial }{\mathit{\theta }}_i}{|}_{{\mathit{\theta }}_i=0}+\frac{\mathit{\partial }{\mathit{\varphi }}_{i,\mathit{\theta }}^s}{\mathit{\partial }{v}_{i,\mathit{\theta }}}(x,z;v_{i,\mathit{\theta }})\frac{\mathit{\partial }{v}_{i,\mathit{\theta }}}{\mathit{\partial }{\mathit{\theta }}_i}{|}_{{\mathit{\theta }}_i=0}.\end{array}$$(22)

From (21(21) $$\begin{array}{c}\hfill {\mathit{\varphi }}_{i,\mathit{\theta }}^s(x,z;v_{i,\mathit{\theta }})=0,(x,z)\in {\mathrm{\Omega }}_{i,\mathit{\theta }},\end{array}$$(21) ), which implies that(23) $$\begin{array}{c}\hfill \frac{\mathit{\partial }{\mathit{\varphi }}_{i,\mathit{\theta }}^s}{\mathit{\partial }{\mathit{\theta }}_i}(x,z;v_{i,\mathit{\theta }}){|}_{{\mathit{\theta }}_i=0}=0,\frac{\mathit{\partial }{\mathit{\varphi }}_{i,\mathit{\theta }}^s}{\mathit{\partial }z}(x,z;v_{i,\mathit{\theta }})\frac{\mathit{\partial }z}{\mathit{\partial }{\mathit{\theta }}_i}{|}_{{\mathit{\theta }}_i=0}=0,\end{array}$$(23)

thus (22(22) $$\begin{array}{c}\hfill \frac{\mathit{\partial }{\mathit{\varphi }}_{i,\mathit{\theta }}^s}{\mathit{\partial }{\mathit{\theta }}_i}(x,z;v_{i,\mathit{\theta }}){|}_{{\mathit{\theta }}_i=0}=\frac{\mathit{\partial }{\mathit{\varphi }}_{i,\mathit{\theta }}^s}{\mathit{\partial }z}(x,z;v_{i,\mathit{\theta }})\frac{\mathit{\partial }z}{\mathit{\partial }{\mathit{\theta }}_i}{|}_{{\mathit{\theta }}_i=0}+\frac{\mathit{\partial }{\mathit{\varphi }}_{i,\mathit{\theta }}^s}{\mathit{\partial }{v}_{i,\mathit{\theta }}}(x,z;v_{i,\mathit{\theta }})\frac{\mathit{\partial }{v}_{i,\mathit{\theta }}}{\mathit{\partial }{\mathit{\theta }}_i}{|}_{{\mathit{\theta }}_i=0}.\end{array}$$(22) ) is simplified to(24) $$\begin{array}{c}\hfill \frac{\mathit{\partial }{\mathit{\varphi }}_{i,\mathit{\theta }}^s}{\mathit{\partial }{v}_{i,\mathit{\theta }}}(x,z;v_{i,\mathit{\theta }}){|}_{{\mathit{\theta }}_i=0}=0.\end{array}$$(24)

Since the operator A is a linear and continuous operator, we have(25) $$\begin{array}{c}\hfill A(\frac{\mathit{\partial }{\mathrm{\Psi }}_{i,\mathit{\theta }}^s}{\mathit{\partial }{v}_{i,\mathit{\theta }}}(x,z;v_{i\mathit{\theta }}){|}_{{\mathit{\theta }}_i=0})=0,\end{array}$$(25)

that is(26) $$\begin{array}{c}\hfill \mathrm{\Delta }(\frac{\mathit{\partial }{\mathrm{\Psi }}_{i,\mathit{\theta }}^s}{\mathit{\partial }{v}_{i,\mathit{\theta }}}{|}_{{\mathit{\theta }}_i=0})+k^2(\frac{\mathit{\partial }{\mathrm{\Psi }}_{i,\mathit{\theta }}^s}{\mathit{\partial }{v}_{i,\mathit{\theta }}}{|}_{{\mathit{\theta }}_i=0})=\mathrm{\Delta }{\mathrm{\Psi }}_i^{\prime }+k^2{\mathrm{\Psi }}_i^{\prime }=0\mathrm{in}\mathrm{\Omega }.\end{array}$$(26)

The Lemma 2.1 is proved. $\square$

Lemma 2.2:

The Fréchet derivative ${\mathrm{\Psi }}_i^{\prime }$ satisfies the following boundary condition:$${\mathrm{\Psi }}_i^{\prime }=-[\frac{\mathit{\partial }({\mathrm{\Psi }}^s+{\mathrm{\Psi }}^{inc})}{\mathit{\partial }z}(\frac{\mathit{\partial }z}{\mathit{\partial }{\mathit{\theta }}_i}{|}_{{\mathit{\theta }}_i=0})]on\mathrm{\Gamma }.$$

Proof:

Solving the derivative of the boundary condition of (17(17) $$\begin{array}{c}\hfill \left\{\begin{array}{c}\mathrm{\Delta }{\mathrm{\Psi }}_{i,\mathit{\theta }}^s+k^2{\mathrm{\Psi }}_{i,\mathit{\theta }}^s=0\mathrm{in}{\mathrm{\Omega }}_{i,\mathit{\theta }},\hfill \\ {\mathrm{\Psi }}_{i,\mathit{\theta }}^s=-{\mathrm{\Psi }}^{inc}\mathrm{on}{\mathrm{\Gamma }}_{i,\mathit{\theta }},\hfill \\ {\mathrm{\Psi }}_{i,\mathit{\theta }}^s(x,z;v_{i,\mathit{\theta }})=\frac{1}{\sqrt{2\mathit{\pi }}}{\int }_{-\infty }^{+\infty }exp(\mathrm{i}x\mathit{\xi }+\mathrm{i}(z-l_{i,\mathit{\theta }})\sqrt{k^2-{\mathit{\xi }}^2})(\mathcal{F}{\mathrm{\Psi }}_{i,\mathit{\theta }Dir}^s)(\mathit{\xi };v_{i,\mathit{\theta }})\mathrm{d}\mathit{\xi }\mathrm{in}{\mathrm{\Omega }}_{l_{i,\mathit{\theta }}},\hfill \end{array}\right.\end{array}$$(17) ) about ${\mathit{\theta }}_i$, we have(27) $$\begin{array}{c}\hfill \frac{\mathit{\partial }{\mathrm{\Psi }}_{i,\mathit{\theta }}^s}{\mathit{\partial }z}(x,z;v_{i,\mathit{\theta }})\frac{\mathit{\partial }z}{\mathit{\partial }{\mathit{\theta }}_i}+\frac{\mathit{\partial }{\mathrm{\Psi }}_{i,\mathit{\theta }}^s}{\mathit{\partial }{v}_{i,\mathit{\theta }}}(x,z;v_{i,\mathit{\theta }})\frac{\mathit{\partial }{v}_{i,\mathit{\theta }}}{\mathit{\partial }{\mathit{\theta }}_i}=-\frac{\mathit{\partial }{\mathrm{\Psi }}^{inc}}{\mathit{\partial }z}(x,z)\frac{\mathit{\partial }z}{\mathit{\partial }{\mathit{\theta }}_i}.\end{array}$$(27)

Hence, from (27(27) $$\begin{array}{c}\hfill \frac{\mathit{\partial }{\mathrm{\Psi }}_{i,\mathit{\theta }}^s}{\mathit{\partial }z}(x,z;v_{i,\mathit{\theta }})\frac{\mathit{\partial }z}{\mathit{\partial }{\mathit{\theta }}_i}+\frac{\mathit{\partial }{\mathrm{\Psi }}_{i,\mathit{\theta }}^s}{\mathit{\partial }{v}_{i,\mathit{\theta }}}(x,z;v_{i,\mathit{\theta }})\frac{\mathit{\partial }{v}_{i,\mathit{\theta }}}{\mathit{\partial }{\mathit{\theta }}_i}=-\frac{\mathit{\partial }{\mathrm{\Psi }}^{inc}}{\mathit{\partial }z}(x,z)\frac{\mathit{\partial }z}{\mathit{\partial }{\mathit{\theta }}_i}.\end{array}$$(27) ), we can get(28) $$\begin{array}{cc}\hfill {\mathrm{\Psi }}_i^{\prime }=\frac{\mathit{\partial }{\mathrm{\Psi }}_{i,\mathit{\theta }}^s}{\mathit{\partial }{v}_{i,\mathit{\theta }}}(x,z;v_{i,\mathit{\theta }})\frac{\mathit{\partial }{v}_{i,\mathit{\theta }}}{\mathit{\partial }{\mathit{\theta }}_i}{|}_{{\mathit{\theta }}_i=0}& =-\frac{\mathit{\partial }{\mathrm{\Psi }}_{i,\mathit{\theta }}^s}{\mathit{\partial }z}(x,z;v_{i,\mathit{\theta }})\frac{\mathit{\partial }z}{\mathit{\partial }{\mathit{\theta }}_i}{|}_{{\mathit{\theta }}_i=0}-\frac{\mathit{\partial }{\mathrm{\Psi }}^{inc}}{\mathit{\partial }z}(x,z)\frac{\mathit{\partial }z}{\mathit{\partial }{\mathit{\theta }}_i}{|}_{{\mathit{\theta }}_i=0}\hfill \\ \hfill & =-[\frac{\mathit{\partial }({\mathrm{\Psi }}^s+{\mathrm{\Psi }}^{inc})}{\mathit{\partial }z}(\frac{\mathit{\partial }z}{\mathit{\partial }{\mathit{\theta }}_i}{|}_{{\mathit{\theta }}_i=0})].\hfill \end{array}$$(28)

The Lemma 2.2 is proved. $\square$

Lemma 2.3:

The Fréchet derivative ${\mathrm{\Psi }}_i^{\prime }$ satisfies the following outgoing radiation condition:$${\mathrm{\Psi }}_i^{\prime }(x,z)=\frac{1}{\sqrt{2\mathit{\pi }}}{\int }_{-\infty }^{+\infty }exp\left(\mathrm{i}x\mathit{\xi }+\mathrm{i}(z-l_H)\sqrt{k^2-{\mathit{\xi }}^2}\right)\left(\mathcal{F}{\mathrm{\Psi }}_{i,Dir}^{{}^{\prime }}\right)(\mathit{\xi })\mathrm{d}\mathit{\xi }in{\mathrm{\Omega }}_{l_H},$$

where $\mathcal{F}{\mathrm{\Psi }}_{i,Dir}^{{}^{\prime }}(\mathit{\xi }):=\mathcal{F}[\frac{\mathit{\partial }{\mathrm{\Psi }}_{i,\mathit{\theta }}^s}{\mathit{\partial }{v}_{i,\mathit{\theta }}}(x,l_H;v_{i,\mathit{\theta }}){|}_{{\mathit{\theta }}_i=0}]\in L^2(\mathbb{R})$, and $l_H$ is chosen such that the half-plane ${\mathrm{\Omega }}_{l_H}\subset {\mathrm{\Omega }}_{i,\mathit{\theta }}$ and ${\mathrm{\Omega }}_{l_H}\subset \mathrm{\Omega }$.

Proof:

This lemma results from the asymptotic behaviour of the scattered field ${\mathrm{\Psi }}_{i,\mathit{\theta }}^s$. Indeed, the scattered field ${\mathrm{\Psi }}_{i,\mathit{\theta }}^s$ satisfies the angular spectrum representation condition as follows:(29) $$\begin{array}{c}\hfill {\mathrm{\Psi }}_{i,\mathit{\theta }}^s(x,z;v_{i,\mathit{\theta }})=\frac{1}{\sqrt{2\mathit{\pi }}}{\int }_{-\infty }^{+\infty }exp\left(\mathrm{i}x\mathit{\xi }+\mathrm{i}(z-l_H)\sqrt{k^2-{\mathit{\xi }}^2}\right)(\mathcal{F}{\mathrm{\Psi }}_{i,\mathit{\theta }Dir}^s)(\mathit{\xi };v_{i,\mathit{\theta }})\mathrm{d}\mathit{\xi }\mathrm{in}{\mathrm{\Omega }}_{l_H},\end{array}$$(29)

where $l_H$ is chosen such that the half-plane ${\mathrm{\Omega }}_{l_H}\subset {\mathrm{\Omega }}_{i,\mathit{\theta }}$ and ${\mathrm{\Omega }}_{l_H}\subset \mathrm{\Omega }$.

Solving the derivative of (29(29) $$\begin{array}{c}\hfill {\mathrm{\Psi }}_{i,\mathit{\theta }}^s(x,z;v_{i,\mathit{\theta }})=\frac{1}{\sqrt{2\mathit{\pi }}}{\int }_{-\infty }^{+\infty }exp\left(\mathrm{i}x\mathit{\xi }+\mathrm{i}(z-l_H)\sqrt{k^2-{\mathit{\xi }}^2}\right)(\mathcal{F}{\mathrm{\Psi }}_{i,\mathit{\theta }Dir}^s)(\mathit{\xi };v_{i,\mathit{\theta }})\mathrm{d}\mathit{\xi }\mathrm{in}{\mathrm{\Omega }}_{l_H},\end{array}$$(29) ) with respect to parameter $v_{i,\mathit{\theta }}$ and taking derivative values at ${\mathit{\theta }}_i=0$, we have(30) $$\begin{array}{cc}& \frac{\mathit{\partial }{\mathrm{\Psi }}_{i,\mathit{\theta }}^s}{\mathit{\partial }{v}_{i,\mathit{\theta }}}{|}_{{\mathit{\theta }}_i=0}\hfill \\ \hfill & =\frac{1}{\sqrt{2\mathit{\pi }}}{\int }_{-\infty }^{+\infty }\{exp\left(\mathrm{i}x\mathit{\xi }+\mathrm{i}(z-l_H)\sqrt{k^2-{\mathit{\xi }}^2}\right)[\frac{\mathit{\partial }}{\mathit{\partial }{v}_{i,\mathit{\theta }}}\left(\mathcal{F}{\mathrm{\Psi }}_{i,\mathit{\theta }Dir}^s\right)(\mathit{\xi };v_{i,\mathit{\theta }})]\}{|}_{{\mathit{\theta }}_i=0}\mathrm{d}\mathit{\xi }\hfill \\ \hfill & =\frac{1}{\sqrt{2\mathit{\pi }}}{\int }_{-\infty }^{+\infty }exp\left(\mathrm{i}x\mathit{\xi }+\mathrm{i}(z-l_H)\sqrt{k^2-{\mathit{\xi }}^2}\right)\mathcal{F}[\frac{\mathit{\partial }{\mathrm{\Psi }}_{i,\mathit{\theta }}^s}{\mathit{\partial }{v}_{i,\mathit{\theta }}}(x,l_H;v_{i,\mathit{\theta }}){|}_{{\mathit{\theta }}_i=0}]\mathrm{d}\mathit{\xi }\hfill \\ \hfill & =\frac{1}{\sqrt{2\mathit{\pi }}}{\int }_{-\infty }^{+\infty }exp\left(\mathrm{i}x\mathit{\xi }+\mathrm{i}(z-l_H)\sqrt{k^2-{\mathit{\xi }}^2}\right)\left(\mathcal{F}{\mathrm{\Psi }}_{i,Dir}^{{}^{\prime }}\right)(\mathit{\xi })\mathrm{d}\mathit{\xi },\mathrm{in}{\mathrm{\Omega }}_{l_H}.\hfill \end{array}$$(30)

The Lemma 2.3 is proved. $\square$

According to the above lemmas, we get the following theorem.

Theorem 2.4:

The Fréchet derivative ${\mathrm{\Psi }}_i^{\prime }$ is the solution of the boundary value problem as follows:(31) $$\begin{array}{c}\hfill \left\{\begin{array}{c}\mathrm{\Delta }{\mathrm{\Psi }}_i^{\prime }+k^2{\mathrm{\Psi }}_i^{\prime }=0in\mathrm{\Omega },\hfill \\ {\mathrm{\Psi }}_i^{\prime }=-[\frac{\mathit{\partial }({\mathrm{\Psi }}^s+{\mathrm{\Psi }}^{inc})}{\mathit{\partial }z}(\frac{\mathit{\partial }z}{\mathit{\partial }\mathit{\theta }}{|}_{\mathit{\theta }=0})],on\mathrm{\Gamma },\hfill \\ {\mathrm{\Psi }}_i^{\prime }(x,z)=\frac{1}{\sqrt{2\mathit{\pi }}}{\int }_{-\infty }^{+\infty }exp(\mathrm{i}x\mathit{\xi }+\mathrm{i}(z-l_H)\sqrt{k^2-{\mathit{\xi }}^2})(\mathcal{F}{\mathrm{\Psi }}_{i,Dir}^{{}^{\prime }})(\mathit{\xi })\mathrm{d}\mathit{\xi },in{\mathrm{\Omega }}_{l_H}.\hfill \end{array}\right.\end{array}$$(31)

where $l_H$ is chosen such that the half-plane ${\mathrm{\Omega }}_{l_H}\subset {\mathrm{\Omega }}_{i,\mathit{\theta }}(i=1,2,\dots ,n)$ and ${\mathrm{\Omega }}_{l_H}\subset \mathrm{\Omega }$.

It is shown by Theorem 2.4 that the Fréchet derivative of operator $\mathcal{K}$ can be characterized by a solution ${\mathrm{\Psi }}^{\prime }$ of the boundary value problem (31(31) $$\begin{array}{c}\hfill \left\{\begin{array}{c}\mathrm{\Delta }{\mathrm{\Psi }}_i^{\prime }+k^2{\mathrm{\Psi }}_i^{\prime }=0in\mathrm{\Omega },\hfill \\ {\mathrm{\Psi }}_i^{\prime }=-[\frac{\mathit{\partial }({\mathrm{\Psi }}^s+{\mathrm{\Psi }}^{inc})}{\mathit{\partial }z}(\frac{\mathit{\partial }z}{\mathit{\partial }\mathit{\theta }}{|}_{\mathit{\theta }=0})],on\mathrm{\Gamma },\hfill \\ {\mathrm{\Psi }}_i^{\prime }(x,z)=\frac{1}{\sqrt{2\mathit{\pi }}}{\int }_{-\infty }^{+\infty }exp(\mathrm{i}x\mathit{\xi }+\mathrm{i}(z-l_H)\sqrt{k^2-{\mathit{\xi }}^2})(\mathcal{F}{\mathrm{\Psi }}_{i,Dir}^{{}^{\prime }})(\mathit{\xi })\mathrm{d}\mathit{\xi },in{\mathrm{\Omega }}_{l_H}.\hfill \end{array}\right.\end{array}$$(31) ).

Theorem 2.5:

The operator $\mathcal{K}$ is Fréchet differentiable and its derivative with respect to parameter $v_i$ is given by the solution of the boundary value problem (31(31) $$\begin{array}{c}\hfill \left\{\begin{array}{c}\mathrm{\Delta }{\mathrm{\Psi }}_i^{\prime }+k^2{\mathrm{\Psi }}_i^{\prime }=0in\mathrm{\Omega },\hfill \\ {\mathrm{\Psi }}_i^{\prime }=-[\frac{\mathit{\partial }({\mathrm{\Psi }}^s+{\mathrm{\Psi }}^{inc})}{\mathit{\partial }z}(\frac{\mathit{\partial }z}{\mathit{\partial }\mathit{\theta }}{|}_{\mathit{\theta }=0})],on\mathrm{\Gamma },\hfill \\ {\mathrm{\Psi }}_i^{\prime }(x,z)=\frac{1}{\sqrt{2\mathit{\pi }}}{\int }_{-\infty }^{+\infty }exp(\mathrm{i}x\mathit{\xi }+\mathrm{i}(z-l_H)\sqrt{k^2-{\mathit{\xi }}^2})(\mathcal{F}{\mathrm{\Psi }}_{i,Dir}^{{}^{\prime }})(\mathit{\xi })\mathrm{d}\mathit{\xi },in{\mathrm{\Omega }}_{l_H}.\hfill \end{array}\right.\end{array}$$(31) ), that is(32) $$\begin{array}{c}\hfill \frac{\mathit{\partial }\mathcal{K}(\overline{v})}{\mathit{\partial }{v}_i}={\mathrm{\Psi }}_i^{\prime }(x,H),i=1,2,\dots ,n.\end{array}$$(32)

4. Multi-parameter identification and shape reconstruction

In this section, we will discuss the inverse scattering problem and reconstruction algorithm. In the inverse problem, given the measured scattering field which is denoted by ${\mathrm{\Psi }}^{test}(x,H)$, we want to determine the unknown parameter vector $\overline{v}$ of fractal rough surface such that(33) $$\begin{array}{c}\hfill \mathcal{K}(\overline{v})-{\mathrm{\Psi }}^{test}(x,H)=0.\end{array}$$(33)

The geometry of the problem is shown in Figure 1 Since $\mathcal{K}$ represents the symbolic relation between the parameters $\overline{v}$ and the scattered fields, it involves coupled integral representations of fields described in Equations (13(13) $$\begin{array}{c}\hfill {\mathrm{\Psi }}^{inc}(\overline{r})+{\int }_{\mathrm{\Gamma }}(\frac{\mathit{\partial }\mathrm{\Phi }(\overline{r},\overline{r^{\prime }})}{\mathit{\partial }{n}^{\prime }}\mathrm{\Psi }(\overline{r^{\prime }})-\mathrm{\Phi }(\overline{r},\overline{r^{\prime }})\frac{\mathit{\partial }\mathrm{\Psi }(\overline{r^{\prime }})}{\mathit{\partial }{n}^{\prime }})\mathrm{d}{s}^{\prime }=\left\{\begin{array}{cc}\mathrm{\Psi }(\overline{r}),\hfill & \overline{r}\in \mathrm{\Omega },\hfill \\ 0,\hfill & \overline{r}\in {\mathbb{R}}^2\backslash \mathrm{\Omega }.\hfill \end{array}\right.\end{array}$$(13) ) and (15(15) $$\begin{array}{c}\hfill {\int }_{-L}^{+L}\mathrm{\Phi }(x,f(x);x^{\prime },f(x^{\prime }))\sqrt{1+{[f^{\prime }(x^{\prime })]}^2}[\frac{\mathit{\partial }\mathrm{\Psi }(\overline{r^{\prime }})}{\mathit{\partial }{n}^{\prime }}{]}_{z^{\prime }=f(x^{\prime })}\mathrm{d}{x}^{\prime }=g^{\mathit{\delta }}(x),x\in [-L,L],\end{array}$$(15) ). The relationship between the rough surface and the scattered filed is non-linear if single scattering is not assumed. In fact, the inversion problem can be written as an optimization problem. That is, find ${\overline{v}}^{\ast }$ such that(34) $$\begin{array}{c}\hfill \Vert \mathcal{K}({\overline{v}}^{\ast })-{\mathrm{\Psi }}^{test}(x,H)\Vert =min\Vert \mathcal{K}(\overline{v})-{\mathrm{\Psi }}^{test}(x,H)\Vert ,\end{array}$$(34)

where $\Vert ·\Vert$ means certain kind of norm.

Figure 1. A schematic of the rough surface inverse scattering problem with the tapered wave incidence.

In the following, a numerical method for the parameters inversion of the rough surface from the measured data is described. For the illumination of the rough surface by M-prescribed incident fields ${\mathrm{\Psi }}_m^{inc}$ ($m=1,2,\dots ,M$), we can obtain the corresponding measurements of scattered fields at nth observation points $(x_n,H)$, and they are denoted by ${\mathrm{\Psi }}_{m,n}^{test}$$(n=1,\dots ,N)$. Based on the measurements ${\mathrm{\Psi }}_{m,n}^{test}$, we wish to evaluate the objective function $C(\overline{v})$ by matching the measurements with the numerical simulations for the trial distributions:(35) $$\begin{array}{c}\hfill C(\overline{v})=(\sum _{n=1}^N|{\mathcal{K}}_{m,n}(\overline{v})-{\mathrm{\Psi }}_{m,n}^{test}{|}^2)/N,m=1,2,\dots ,M,\end{array}$$(35)

where ${\mathcal{K}}_{m,n}(\overline{v})={\mathrm{\Psi }}_m^s(x_n,H)$ represent the scattered fields which are dependent on $\overline{v}$. The dimension of vector $\overline{v}$ is determined by the assumed mathematical model of the surface profile.

To solve the optimization problem (35(35) $$\begin{array}{c}\hfill C(\overline{v})=(\sum _{n=1}^N|{\mathcal{K}}_{m,n}(\overline{v})-{\mathrm{\Psi }}_{m,n}^{test}{|}^2)/N,m=1,2,\dots ,M,\end{array}$$(35) ) for the parameters inversion of the rough surface function f by BFGS method, we need to consider the existence of the gradient $\mathrm{\nabla }C(\overline{v})$ of the objective function $C(\overline{v})$. From (35(35) $$\begin{array}{c}\hfill C(\overline{v})=(\sum _{n=1}^N|{\mathcal{K}}_{m,n}(\overline{v})-{\mathrm{\Psi }}_{m,n}^{test}{|}^2)/N,m=1,2,\dots ,M,\end{array}$$(35) ) and Theorem 2.5, the following theorem can be proved.

Theorem 3.1:

The gradient $\mathrm{\nabla }C(\overline{v})$ of the objective function $C(\overline{v})$ exists.

Based on the above analysis, the reconstruction algorithm is summarized as follows:

• Step 1. Choose an initial guess of ${\overline{v}}_0$, an accuracy $\mathit{\epsilon }>0$ and a symmetric positive definite matrix $B_0$. Set $k=0$.

• Step 2. Solve the boundary value problem (31(31) $$\begin{array}{c}\hfill \left\{\begin{array}{c}\mathrm{\Delta }{\mathrm{\Psi }}_i^{\prime }+k^2{\mathrm{\Psi }}_i^{\prime }=0in\mathrm{\Omega },\hfill \\ {\mathrm{\Psi }}_i^{\prime }=-[\frac{\mathit{\partial }({\mathrm{\Psi }}^s+{\mathrm{\Psi }}^{inc})}{\mathit{\partial }z}(\frac{\mathit{\partial }z}{\mathit{\partial }\mathit{\theta }}{|}_{\mathit{\theta }=0})],on\mathrm{\Gamma },\hfill \\ {\mathrm{\Psi }}_i^{\prime }(x,z)=\frac{1}{\sqrt{2\mathit{\pi }}}{\int }_{-\infty }^{+\infty }exp(\mathrm{i}x\mathit{\xi }+\mathrm{i}(z-l_H)\sqrt{k^2-{\mathit{\xi }}^2})(\mathcal{F}{\mathrm{\Psi }}_{i,Dir}^{{}^{\prime }})(\mathit{\xi })\mathrm{d}\mathit{\xi },in{\mathrm{\Omega }}_{l_H}.\hfill \end{array}\right.\end{array}$$(31) ) by the integral equation method to get the numerical solutions $\mathrm{\nabla }C({\overline{v}}_k)$. If $\Vert \mathrm{\nabla }C({\overline{v}}_k)\Vert <\mathit{\epsilon }$, stop the algorithm. Otherwise, compute the direction ${\overline{p}}_k$ by$${\overline{p}}_k=-B_k^{-1}\mathrm{\nabla }C({\overline{v}}_k).$$

• Step 3. Find an acceptable ${\mathit{\alpha }}_k$ such that$$C({\overline{v}}_k+{\mathit{\alpha }}_k{\overline{p}}_k)=\underset{\mathit{\alpha }\ge 0}{min}C({\overline{v}}_k+\mathit{\alpha }{\overline{p}}_k).$$

• Step 4. Set ${\overline{v}}_{k+1}={\overline{v}}_k+{\mathit{\alpha }}_k{\overline{p}}_k$, then compute$$\begin{array}{cc}\hfill {\overline{s}}_k& ={\mathit{\alpha }}_k{\overline{p}}_k,\hfill \\ \hfill {\overline{y}}_k& =\mathrm{\nabla }C({\overline{v}}_{k+1})-\mathrm{\nabla }C({\overline{v}}_k),\hfill \\ \hfill B_{k+1}& =B_k+\frac{{\overline{y}}_k{\overline{y}}_k^{\mathrm{T}}}{{\overline{y}}_k^{\mathrm{T}}{\overline{s}}_k}-\frac{B_k{\overline{s}}_k{\overline{s}}_k^{\mathrm{T}}{B}_k}{{\overline{s}}_k^{\mathrm{T}}{B}_k{\overline{s}}_k}.\hfill \end{array}$$

• Step 5. Let $k=k+1$, and repeat Steps 2–4.

In general, the objective function is a non-linear function of the parameters. It is difficult to predict the behaviour of the objective functions with respect to variations in these parameters. The BFGS method builds up the solution of the optimization problem by successively solving the direct scattering problem (to get $C(\overline{v})$) and the adjoint problem (to get $\mathrm{\nabla }C(\overline{v})$). The proposed algorithm generates a sequence of values $\{{\overline{v}}_k\}$ approximate to the minimum ${\overline{v}}^{\ast }$ of the objective function $C(\overline{v})$. This implies that for each iteration, several runs of the fast forward solver are required to retrieve the surface profile in realistic time, either to compute the objective function or to compute the gradient of the objective function.

5. Numerical results

In this section, we illustrate the performance of the inversion algorithm. In the following numerical examples, we consider two cases of reconstruction, respectively. One is the reconstruction with a single incident wave; the other one is the reconstruction with multi-angle and multi-frequency incident waves.

5.1. Single incidence reconstruction

We consider the inverse scattering problem for the fractal rough surface $\mathrm{\Gamma }$ which is generated by the surface height function f(x), where the parameters $M_f=9,b=2.71828,{\mathit{\phi }}_n=0$, and the other parameters $\widehat{h},\widehat{D}$, K are given for each numerical experiment below. The surface is truncated in a finite integration interval $[-L,L]$$(L=80)$. For tapered wave incidence, we set the wavelength $\mathit{\lambda }$ and the incidence angle ${\mathit{\theta }}_{inc}$. To avoid the edge effect and assure the computational precision, we set the tapered factor g = L / 2 as in [29]. For solving the direct problem, the rough surface is divided into N segments ($N=\frac{2L}{\mathrm{\Delta }x}$), where $\mathrm{\Delta }x=\mathit{\lambda }/10$ is the space interval. The scattered field measurements ${\mathrm{\Psi }}_{m,n}^{test}$ are obtained at the angles ${\mathit{\theta }}_{s,k}=-\frac{\mathit{\pi }}{4}+\frac{k\mathit{\pi }}{60}$, $(k=0,1,\dots ,30)$ and $H=10$. They are applied to replace the real collected GPR-type data to test the reconstruction algorithm. Let ${\widehat{h}}_{num}$, ${\widehat{D}}_{num}$ and $K_{num}$ denote the inversion results of $\widehat{h},\widehat{D},K$ respectively, $N_{it}$ denote the number of iterations, $Ma{x}_{er}$ denote the maximum error between the reconstructed surface $\widehat{f}(x)$ and the original real one f(x), i.e. $Ma{x}_{er}:={max}_{x\in [-L,L]}|f(x)-\widehat{f}(x)|$.

For two parameters inversion, set ${\mathit{\theta }}_{inc}=\frac{\mathit{\pi }}{6},\mathit{\lambda }=1$ and $K=0.3$. We show the inversion results of $\widehat{h}$ and $\widehat{D}$ under different conditions in Tables 1–5. As we can see clearly from the Table 1, although for these examples we all choose the initial data as ${\widehat{h}}_0=0.2,{\widehat{D}}_0=1.2$, our approach can reconstruct the parameters of the rough surface very well.

Figure 2. Reconstruct the surface for $\widehat{h}=0.04,\widehat{D}=1.3$ (left) and $\widehat{h}=0.1,\widehat{D}=1.3$ (right) with 1% noise.

Figure 3. Reconstruct the surface for $(\widehat{h})=0.3$, $(\widehat{D})=1.3$ (left) and $(\widehat{h})=0.4$, $(\widehat{D})=1.5$ (right) with 1% noise.

Table 1. The two-parameter $(\widehat{h},\widehat{D})$ inversion results of the fractal rough surface.

$\widehat{h}$	$\widehat{D}$	${\widehat{h}}_{num}$	${\widehat{D}}_{num}$	$N_{it}$	$Ma{x}_{er}$
0.04	1.3	0.0401	1.3059	26	$<{10}^{-3}$
0.1	1.3	0.0999	1.2981	14	$<{10}^{-3}$
0.3	1.3	0.3001	1.3003	4	$<{10}^{-3}$
0.4	1.5	0.4003	1.5001	21	$<{10}^{-3}$

In fact, there are errors for the measurements. To test the influence of measurement error, some random noises are added on the measurement information. For the relative noise level $\mathit{\delta }$, the measurements are constructed in the following way:(36) $$\begin{array}{c}\hfill {\mathrm{\Psi }}_{\mathit{\delta }}^{test}={\mathrm{\Psi }}^{test}+\mathit{\delta }(\frac{\mathit{\vartheta }}{{|\mathit{\vartheta }|}_{\infty }})\circ {\mathrm{\Psi }}^{test},\end{array}$$(36)

where $\mathit{\vartheta }$ is a random vector, the element of the vector ${\mathit{\vartheta }}_i\in$$\mathbf{N}(0,1)$ normally distributed, ${|\mathit{\vartheta }|}_{\infty }={max}_i|{\mathit{\vartheta }}_i|$. For two Matrices, A and B, of the same dimensions, the Hadamard product $A\circ B$ is a matrix of the same dimensions, and $A\circ B=[a_{i,j}{b}_{i,j}]$.

So we show the inversion results with the measurement data which bring 1% noise for different parameters of the rough surface in Table 2. We can see that the numerical results are close to the theoretical values. Meanwhile, the results of shape reconstruction of the fractal rough surfaces are also given in Figures 2 and 3.

Table 2. The two-parameter $(\widehat{h},\widehat{D})$ inversion results of the fractal rough surface with 1% noise level.

Noise (%)	$\widehat{h}$	$\widehat{D}$	${\widehat{h}}_{num}$	${\widehat{D}}_{num}$	$N_{it}$	$Ma{x}_{er}$
1	0.04	1.3	0.0402	1.2949	12	$<{10}^{-3}$
1	0.1	1.3	0.1000	1.2944	28	$<{10}^{-3}$
1	0.3	1.3	0.2999	1.2996	11	$<{10}^{-3}$
1	0.4	1.5	0.4005	1.5003	16	$<{10}^{-3}$

In Tables 3–5, we analyse the influence of some parameters for the algorithm as mentioned above, such as the error of the measurements, the initial guesses and the survey positions. The comparisons of the different errors of the measurement data for the numerical results are shown in Table 3. The maximum error is less than ${10}^{-3}$ when the noise level does not exceed 5%. Based on the above comparisons, we can see that the results are satisfactory.

Table 3. The two-parameter $(\widehat{h},\widehat{D})$ inversion results of the fractal rough surface with 1–5% noise level.

Noise (%)	$\widehat{h}$	$\widehat{D}$	${\widehat{h}}_{num}$	${\widehat{D}}_{num}$	$N_{it}$	$Ma{x}_{er}$
1	0.3	1.3	0.2999	1.2996	11	$<{10}^{-3}$
2	0.3	1.3	0.3002	1.2986	5	$<{10}^{-3}$
3	0.3	1.3	0.3001	1.2993	9	$<{10}^{-3}$
4	0.3	1.3	0.3000	1.3012	9	$<{10}^{-3}$
5	0.3	1.3	0.3000	1.2986	8	$<{10}^{-3}$

In Table 4, we can see the influence of initial guesses on the numerical results. It is found that when the initial guesses are not far from the real one, the inversion results are in fairly good agreement with the theoretical values of these numerical examples. In this paper, we choose some initial parameters to test the algorithm. In practical problems, the better initial guesses could be obtained by some analytical approximation methods or some direct inversion methods (see, e.g. [16,17]). From Table 5, we can see that the survey position has little influence on the inversion results. The calculation results are in fairly good agreement with different survey positions for all these numerical examples.

These parameters of the rough surface have a very complex impact on the scattering field, which together determine the distribution of scattered field. From the inversion results which are shown in the Tables 1–5, we can see that some information of the rough surface, such as parameters $\widehat{h},\widehat{D}$, can be retrieved from the measurements very well, and the parameters such as the error of the measurements, the initial data and the survey position have different effects on the simulation results. The numerical results show that the algorithm yields satisfactory reconstructions for some rough surface profiles.

Table 4. The two-parameter $(\widehat{h},\widehat{D})$ inversion results of the fractal rough surface with different initial values.

$({\widehat{h}}_0,{\widehat{D}}_0)$	$\widehat{h}$	$\widehat{D}$	${\widehat{h}}_{num}$	${\widehat{D}}_{num}$	$N_{it}$	$Ma{x}_{er}$
(0.2, 1.2)	0.3	1.3	0.3001	1.3003	4	$<{10}^{-3}$
(0.1, 1.2)	0.3	1.3	0.3001	1.3017	8	$<{10}^{-3}$
(0.1, 1.1)	0.3	1.3	0.3004	1.2998	7	$<{10}^{-3}$
(0.2, 1.8)	0.3	1.3	0.2998	1.3006	10	$<{10}^{-3}$
(0.4, 1.5)	0.3	1.3	0.3003	1.3009	6	$<{10}^{-3}$

Table 5. The two-parameter $(\widehat{h},\widehat{D})$ inversion results of the fractal rough surface with different survey positions.

H	$\widehat{h}$	$\widehat{D}$	${\widehat{h}}_{num}$	${\widehat{D}}_{num}$	$N_{it}$	$Ma{x}_{er}$
10	0.3	1.3	0.3001	1.3003	4	$<{10}^{-3}$
1	0.3	1.3	0.2998	1.2998	10	$<{10}^{-3}$
100	0.3	1.3	0.2996	1.3010	5	$<{10}^{-3}$

For three parameters inversion, we still let ${\mathit{\theta }}_{inc}=\frac{\mathit{\pi }}{6},\mathit{\lambda }=1$, and the initial guess ${\widehat{h}}_0=0.2,{\widehat{D}}_0=1.2,K_0=0.4$, $B_0=diag(6,1,6)$. Note that the estimated parameters are used as initial guesses to start the algorithm again if the error value Max${}_{er}>{10}^{-3}$ after more than 30 iterations. In Table 6, we found that despite the relatively large number of iterations used, the approach can reconstruct the parameters of the rough surface as well. The surface profile of the fractal rough surface is also reconstructed in Figure 4.

Table 6. The three-parameter $(\widehat{h},\widehat{D},K)$ inversion results of the fractal rough surface.

$\widehat{h}$	$\widehat{D}$	K	${\widehat{h}}_{num}$	${\widehat{D}}_{num}$	$K_{num}$	$N_{it}$	$Ma{x}_{er}$
0.1	1.3	0.5	0.0999	1.2904	0.5000	$30\ast 2+20$	$<{10}^{-3}$

Figure 4. Reconstruct the surface for three parameters $\widehat{h}=0.1,\widehat{D}=1.3,K=0.5.$

5.2. Multi-angle and multi-frequency reconstruction

We now move to consider the incidence of multi-angle and multi-frequency. Firstly, we consider the multiple incidence strategy with a single frequency. For one incident angle, we use the scattered field measurements at all receivers simultaneously and then change the objective function after each determination of the parameters of rough surfaces corresponding to each incident angle. Secondly, we consider the multi-frequency strategy with a single incident angle. Similar to the multiple incidence strategy with the single frequency, for one fixed frequency, we use the scattered field measurements at all receivers simultaneously and then change the objective function after each determination of the parameters of rough surfaces corresponding to each frequency. To improve the results of reconstruction, the combined approach can be implemented as the multi-frequency strategy followed by multiple incidence strategy (which is implemented by the multiple incidence strategy for each frequency from low to high). Conversely, it also can be implemented as the multiple incidence strategy followed by multi- frequency strategy (which is implemented the multiple frequency strategy for each incident angle from small to big).

Figure 5. Reconstruct the surface for $\widehat{h}=0.1,\widehat{D}=1.5$ with 1% noise by single-incidence (The top), multiple-angle incidence (The middle) and multiple-frequency incidence (The bottom).

Figure 6. Reconstruct the surface for three parameters $\widehat{h}=0.1,\widehat{D}=1.3,K=0.5.$ with 1% noise.

Let the incident angles ${\mathit{\theta }}_{inc}$ = $\frac{\mathit{\pi }}{6}$, $\frac{\mathit{\pi }}{4}$ and $\frac{\mathit{\pi }}{3}$, respectively, and set $\mathit{\lambda }=1$, $K=0.3$. At ${\mathit{\theta }}_{inc}$ = $\frac{\mathit{\pi }}{6}$, the inversion algorithm starts with the initial guess ${\widehat{h}}_0=0.2,{\widehat{D}}_0=1.2$ for certain inversion iterations. The above inversion results are then used as initial values in the inversion algorithm at ${\mathit{\theta }}_{inc}=\frac{\mathit{\pi }}{4}$ for certain inversion iterations. This process is repeated at the incident angle ${\mathit{\theta }}_{inc}=\frac{\mathit{\pi }}{3}$. The numerical inversion results are listed in the Table 7 to show a comparison between the inversion results with multiple incidence (M-A) and single incidence (S-I) at a fixed frequency ($\mathit{\lambda }=1$), where the incident angle of S-I is $\frac{\mathit{\pi }}{6}$. All these results are obtained for the same objective function. The multiple incidence strategy need only five inversion iterations at each incident angle. It can achieve the same precision with fewer iterations. We can see that the multiple incidence strategy shows a better result in the parameters reconstruction.

Table 7. The two-parameter $(\widehat{h},\widehat{D})$ inversion results with multiple incidence and single incidence at a fixed frequency.

Inc	Noise (%)	$\widehat{h}$	$\widehat{D}$	${\widehat{h}}_{num}$	${\widehat{D}}_{num}$	$N_{it}$	$Ma{x}_{er}$
M-A	0	0.1	1.5	0.1000	1.5000	$5\ast 3$	$<{10}^{-6}$
S -I	0	0.1	1.5	0.1000	1.4995	35	$<{10}^{-6}$
M-A	1	0.1	1.5	0.1000	1.5020	$5\ast 3$	$<{10}^{-3}$
S -I	1	0.1	1.5	0.0987	1.4421	100	$>{10}^{-3}$

In Table 8, we show the parameters inversion results as $K=0.3$ with different wavelengths $\mathit{\lambda }=4,1$ and 0.5, respectively, for incident angle ${\mathit{\theta }}_{inc}=\frac{\mathit{\pi }}{6}$. We start the algorithm at $\mathit{\lambda }=4$ with the initial guess of the parameters ${\widehat{h}}_0=0.2,{\widehat{D}}_0=1.2$ for certain inversion iterations. The inversion results are then used as the initial values in the inversion algorithm at $\mathit{\lambda }=1$ for certain inversion iterations. This process is repeated at the wavelength $\mathit{\lambda }=0.5$, where M-F represent the multi-frequency incidence. By the reconstructing parameters $h_{num}$ and $D_{num}$, we also show the inversion results of the rough surface for different cases in Figure 5. Based on the numerical results, we can see that the results of reconstruction can be improved with multi-frequency strategy. It requires fewer iterations to achieve the same precision.

Table 8. The two-parameter $(\widehat{h},\widehat{D})$ inversion results of the fractal rough surface with multiple frequency incidence.

Inc	Noise (%)	$\widehat{h}$	$\widehat{D}$	${\widehat{h}}_{num}$	${\widehat{D}}_{num}$	$N_{it}$	$Ma{x}_{er}$
M-F	0	0.1	1.5	0.1000	1.5000	$15\ast 2$	$<{10}^{-6}$
M-F	1	0.1	1.5	0.0999	1.5021	$15\ast 3$	$<{10}^{-3}$

Table 9. The three-parameter $(\widehat{h},\widehat{D},K)$ inversion results of the fractal rough surface.

$\widehat{h}$	$\widehat{D}$	K	Noise	${\widehat{h}}_{num}$	${\widehat{D}}_{num}$	$K_{num}$	$N_{it}$	$Ma{x}_{er}$
0.1	1.3	0.5	1%	0.1002	1.2966	0.5000	$30+27$	$<{10}^{-3}$

For three parameters inversion, multiple incident angles and multiple frequency can be used simultaneously. The initial guesses ${\widehat{h}}_0=0.2,{\widehat{D}}_0=1.2$, $K_0=0.4$ and $B_0=diag(6,1,6)$ are used to start the algorithm. Note that the estimated parameters are used as initial guess to start the algorithm again for other incident angles or frequencies if the error value Max${}_{er}>{10}^{-3}$ after more than 30 iterations. In Table 9, we only used one wave length $(\mathit{\lambda }=4)$ and two incident angles $({\mathit{\theta }}_{inc}=\frac{\mathit{\pi }}{6},\frac{\mathit{\pi }}{4})$ simultaneously; we found that despite the relatively large number of iterations used, the approach can reconstruct the parameters of the rough surface as well. The surface profile of the fractal rough surface is also given in Figure 6.

6. Summary and conclusions

We present an inversion algorithm to reconstruct the parameters of an unbounded fractal rough surfaces with tapered wave incidence. Given the scattered field measurements, we reconstruct the parameters of unbounded rough surfaces. First, the parameters inversion problem is viewed as an optimization problem, and the unknown parameters of the rough surfaces are estimated by minimizing an objective function. By the related Fréchet derivative with respect to the surface parameters of the solution of direct scattering problems, an efficient reconstruction algorithm is proposed. Second, the influence of measurement error, initial data, survey position and different rough surfaces on the inversion results is studied. Finally, multiple incidence strategy is used to improve the conditioning of the reconstruction problem. The incidence strategy of multi-angle and multi-frequency can make us obtain better reconstruction results than single incidence with much less iterations. The numerical results show that the proposed algorithm is accurate, and the inversion results can be obtained with reasonable computing times.

Acknowledgements

The authors would like to thank the referees for their constructive comments and suggestions which helped improve the paper.

Additional information

Funding

This work was partly supported by the National Natural Science Foundation of China [grant number 11401183], [grant number 11271113]; China Postdoctoral Science Foundation funded project [grant number 2015M571858]; and Natural Science Foundation of Heilongjiang Province of China [grant number A2015008].

Notes

No potential conflict of interest was reported by the authors.