Reconstruction of rough surfaces from a single receiver at grazing angle

The paper develops a method for recovering a one-dimensional rough surface profile from scattered wave field, using a single receiver and repeated measurements when the surface is moving with respect to source and receiver. This extends a previously introduced marching method utilizing low grazing angles, and addresses the key issue of the requirement for many simultaneous receivers. The algorithm recovers the surface height below the receiver point step-by-step as the surface is moved, using the parabolic wave integral equations. Numerical examples of reconstructed surfaces demonstrate that the method is robust in both Dirichlet and Neumann boundary conditions, and with respect to different roughness characteristics and to some degree of measurement noise.


Introduction
The calculation of wave scattering by rough surfaces arises in numerous physical problems [1][2][3][4][5][6], and the inverse problem of recovering rough surfaces is a crucial goal in a wide range of engineering applications, ranging from geophysical remote sensing and ocean acoustics (including under-ice monitoring) to nondestructive testing [7,8].Many works have focused on rough surface reconstruction via scattered data, e.g. .Such approaches include Kirchhoff [25], Rytov approximation [9], iterative [19,21] and in particular Newton methods [15][16][17][18], linear sampling [13,14], and more recently deep learning [24].Most of these use single frequency data while others have utilized multiple frequencies [11].In particular, in [29][30][31][32][33], efficient reconstruction algorithms have been developed under the regime of grazing angle incidence using the parabolic wave equation.When the incident angle is small, and forward propagation dominates, the Helmholtz equation can be approximated by the parabolic wave equation.This allows the boundary integral equations to be replaced by the Parabolic Integral Equations (PIE) [34], using which the computational cost of both direct and inverse problems can be significantly reduced.However a drawback of some of the above methods is that in practical applications the scope for number of receivers may be limited, and indeed their placement may interfere with the scattered field to be measured.
One approach to this issue is to allow the measurement apparatus to move with respect to the surface (or equivalently fixed source and receiver with a moving surface) and take repeated measurements using fewer receivers.Such a situation is feasible for applications such as NDT of surface roughness, or under-ice monitoring.However, this poses additional complexity for the computational implementation.
In this paper, an algorithm is developed to reconstruct the rough surface which is moving, for scattering at the grazing angle.The algorithm relies on the scattered data at a single fixed point, thus only one receiver is needed for the reconstruction.This approach is analogous to previous works [29,30] based on PIE, but in those cases, N receivers were required where N is number of unknowns determined by the discretisation size.Both the Dirichlet (TE incident field) and Neumann (TM) boundary conditions are treated here.In this marching approach, the surface profile is recovered point by point as the surface moves.We first recover the surface current with which the surface height below the receiver point can then be reconstructed using the coupled integral equations.
A variety of numerical examples are presented to validate the proposed method.The source is taken by the Gaussian beam with a single frequency.Two types of rough surfaces are tested: Gaussian and sub-fractal.It is found that the proposed method can successfully capture the detailed features of the rough surface.The method is tested with noisy data; although this introduces fine-scale oscillations, reconstruction still agrees well with the actual surface, and these oscillations can easily be filtered out.The error between the reconstructed and actual surfaces in the ℓ 2 -norm is measured with respect to the discretisation size, surface type, surface scale, and the receiver height.
The layout of the paper is structured as follows.In section 2, we give a brief overview of the direct problem using parabolic integral equations.The problem setting and reconstruction algorithm are described in section 3, which includes the detailed formulations for both the Dirichlet and Neumann boundary conditions.Numerical tests are carried out in section 4, where the method is examined by both the error and the performance with noisy data.Finally, concluding remarks and future work are discussed in section 5.

Mathematical Background
Boundary integral formulations of rough surface scattering problem have been extensively studied, for the full Green's function (Helmholtz regime), for example in [2,3], and for the parabolic integral equations (PIE) on which the current study is based [34][35][36].We briefly review the assumptions and the governing integral equations for the parabolic regime studied here.Consider a scattering problem in which an electromagnetic field is incident at a low grazing angle on a perfectly conducting corrugated surface varying in one dimension only, so that the problem becomes two-dimensional and scalar.The incident field is taken to be time-harmonic and either horizontally or vertically plane polarized, i.e. transverse electric (TE) or transverse magnetic (TM).Time dependence may be suppressed so that we consider the time-reduced component.The coordinate axes are x and z, where x is horizontal, and z is vertical.We denote the time-harmonic wave field as E(x, z).The governing equation for the wave field is the Helmholtz equation, where k is the wavenumber.If the incident angle is small enough so that the forward propagation dominates, then there is a slowly varying component of the wave field, denoted by ψ(x, z), given by ψ(x, z) = E(x, z) exp(−ikx). ( With this approximation, the backscattered component is neglected, and the slowly varying wave field can be governed by the parabolic wave equation: The Green's function for the parabolic wave equation is in the form of when x < x, and G = 0 otherwise [34,36].We denote the incident wave field, the scattered wave field, and the total wave field as ψ i (x, z), ψ s (x, z), and ψ(x, z), respectively, so that The forward problem (to obtain the scattered field from a known surface) is typically solved via the Kirchhoff-Helmholtz equations which relate the scattered and incident fields via the surface currents (J or K depending on polarization).The inverse problem tackled below is the reconstruction of a rough surface from scattered data using a single receiver, by allowing the surface to move with respect to the source and receiver.Consider first an arbitrary rough surface z = h(x) on the domain [0, L].The source will be a Gaussian beam propagating at a relatively small angle to the horizontal, so we can consider the slowly varying component, and replace the Kirchhoff-Helmholtz equation by the coupled parabolic integral equations relating the incident field and the scattered field [36].Let the space V contain all the points lying above the surface with V = {(x, z) : z > h(x)}, and let the space S contain the points lying on the surface with S = {(x, z) : z = h(x)}.Two boundary conditions are considered in this paper, Dirichlet and Neumann.If we assume the Dirichlet boundary condition, which corresponds to a transverse electric (TE) field impinging on a perfectly conducting surface, then the coupled integral equations are where points r and r both lie on the surface, and where r again lies on the surface, and r is any point above the surface.For simplicity, we denote ψ ′ := ∂ψ/∂z.On the other hand, if we assume the Neumann boundary condition, which corresponds to a transverse magnetic (TM) polarized field or an acoustically hard surface, then the coupled integral equations become where points r and r both lie on the surface, and where the point r still lies on the surface, but point r can be any point in the medium.

S R Reconstruction Domain Moving Direction
Figure 1: Problem setting for reconstruction of rough surface moving with respect to source and receiver, in which the source (S) and the receiver (R) are located at x = 0 and x = L, respectively.The reconstruction domain is the region between the source and receiver planes.
3 Reconstruction of moving surface

Problem setting
This paper deals with the situation where the surface moves horizontally in the negative x direction with respect to a fixed source and receiver (see fig. 1).This is obviously equivalent to source and receiver moving over a static surface, and the formulation here is chosen for convenience.
The rough segment to be recovered is considered as part of an extended surface with a flat part and rough part.For simplicity the length of each of these regions is set to be L, making the total length of the surface 2L.The source and receiver are above the surface, at given heights and separated horizontally by distance L. Define the reconstruction domain to be the region lying between the vertical planes of source and receiver, x = 0, L respectively.An illustrative figure of the physical setting is shown in Fig. 1.
Thus the surface profile h(x) is initially flat for x ≤ L and irregular for x ≥ L. The region [0, 2L] is discretized using 2N + 1 equally spaced points x n = n∆x for n = 0, 1, • • • , 2N where ∆x = L/N .The surface values at initial time in the area [0, 2L] can be thought of the discretized values in which negative q corresponds to the flat part, and positive q corresponds to the rough part of the extended surface.There are in total N + 1 time steps t 0 to t N .In the inverse problem from time t j to t j+1 , the surface moves by one spatial distance ∆x towards the left and a new measurement of scattered field is taken.
The algorithm is applied within the reconstruction domain [0, L] ≡ [x 0 , x N ] between the source and receiver.As surface moves, the surface height depends on the time step.we denote the surface at the initial time step t 0 as z = h 0 (x) with x ∈ [0.2L].The repositioned surfaces in the reconstruction domain [0, L] at the time step t j for j ≥ 1 is denoted by h j (x) given by Accordingly, the spaces V (t j ) and S(t j ) represent the points lying above and on the surface, respectively, at the step t j , which are given by V (t j ) = {(x, z) : z > h j (x), x ∈ [0, L]}, and

S(t
It is clear that at step t j , the last point of surface in the reconstruction domain is h j (L) = h j (x N ) ≡ H j .Finally, we denote by ψ ′ j , ψ j , and ψ s j the surface currents and scattered field at time step t j .

Reconstruction in the Dirichlet case
Assume first that the Dirichlet boundary condition holds.For the rough surface at the reconstruction domain at step t j , eq. ( 5) and eq.( 6) become and Returning to the discretized problem, surface values H n for n ∈ {−N, .., 0} correspond to the flat segment which lies initially within the reconstruction interval [0, L].The algorithm seeks to recover the unknown surface values H n for n ∈ {1, ...N }.At step t 1 all surface values within the interval are known except that at x = L which is now H 1 ≡ h 1 (x N ).This value is recovered as described below, and the procedure is carried out repeatedly.More generally, at each step t j , only the right-most value H j of those on the interval [0, L] is unknown and is to be recovered.At each step, the reconstruction algorithm is an extension of the marching procedure in the previous works [29,31] as we now describe.The method requires a sequence of scattered field measurements ψ s j (r At time step t j , the surface values H 0 , • • • , H j−1 are known.For any point r n = (x n , h j (x n )) ∈ S(t j ) with x n ∈ [x 1 , x N ], eq. ( 10) is written as a sum of n subintegrals with We assume that the surface derivative ψ ′ , and the Green's function (except in the last summand where it is singular) vary slowly over the subintervals, and can be treated as constant.For n > 1 the integral formula then becomes At this point, we assume that the surface value at the last point is zero, By assumption all surface values required to evaluate the sum are known.If we take n = 1, 2, • • • , N , a linear system can be generated with where A D j is a matrix of size N × N given by for 1 ≤ m ≤ n ≤ N , and Ψ ′ j and Ψ i j are vectors of size N with By solving the linear system, the surface derivative on all the discrete surface points over the reconstruction domain can be obtained.The singularity in the diagonal term of A D j has been treated in [29,Appendix B], which is given by A where γ = k/2|h ′ j (x n )|, and C and S are the cosine and sine Fresnel integrals.Note that the system A D j is lower-triangular, thus the inversion of this system is computationally cheap.
We use the second integral equation (eq.( 6)) to recover the surface height at the end h j (x N ).The receiver is at the point where ζ is a fixed value above the rough surface.Under the assumption that the surface derivative and Green's function can be viewed as constant over the subinterval, the scattered field at the point of receiver at time step t j can be written by For the last integral, we approximate the integral using the midpoint: where ξ N = 1/2(x N −1 + x N ), and α = 1/2(i/2πk) 1/2 .As we know the surface derivative ψ ′ at all the discrete points, and all the surface values other than the last one h(x N ), then the first (N − 1) terms in eq. ( 17) can be calculated.To recover the surface height h j (x N ), the problem becomes solving an exponential equation.The logarithm of a complex-valued number is ambiguous, for z ∈ C, log z = ln(|z|)+iarg(z), where arg(z) = Arg(z)+2kπ with Arg(z) ∈ [−π, π] and k ∈ Z.However, if we assume that the surface height is small enough, we can simply take the principal component, namely, we use Ln(z) = ln|z| + iArg(z).
Combining the two equations eq. ( 17) and eq.( 18), the surface height at the last point is recovered with where Finally, at time step t j , the last point surface height is equivalent to H j , and we set H j := h j (x N ).

Reconstruction in the Neumann case
Similar treatment can be applied to the Neumann case.We still first assume that H n = h 0 (x N ) := 0. Under the assumption that the total surface wave and the derivative of the Green's function vary slowly over each subinterval and can be treated as constants, then at any time step t j , eq. ( 7) can be written by N subintegrals for the points r n = (x n , h j (x n )), r m = (x m , h j (x m )), and r = (x, h j (x)).A linear system is then generated with where A N j is a matrix of size N × N given by for 1 ≤ m ≤ n ≤ N , and Ψ j is the vector of total wave along surface of size N given by and Ψ i j is the vector of incident wave along surface given in eq. ( 16).The singular integral in the diagonal term of A N j has been calculated in [31], with where β = −i/2 ik/2π.The surface current, Ψ j , can be obtained by inverting the lowtriangular matrix A N j .We use eq.( 8) to reconstruct the surface profile to the right, which is h(x N ).At the point of receiver r ζ N = (x N , ζ), eq. ( 8) can be expressed in the form: where the terms r m , r N , r are again evaluated on the surface S(t j ).As in the Dirichlet case, the first (N − 1)-sum is known.We directly approximate the last integral with where ξ N is the midpoint defined by ξ N = (x N + x N −1 )/2.Rearrange eq. ( 25), and take the absolute value both sides, we have Provided that the receiver is higher than the surface, the surface height h j (X N ) can be obtained via where Finally, at the time step t j , we set H j := h j (x N ).

Numerical results
The performance of the proposed method is examined in this section.For ease of comparison with previous PE results, in the numerical examples the wave number is set to be k = 1, resulting in the wavelength of λ = 2π.The extended surface is made up of a flat surface and a rough surface, which in total have length 2L, where L = 200, approximately 31.8λ.The source is a Gaussian beam, given by where it is centred at (x 0 , z 0 ) with z 0 = 22.4,and the beam width w is given by w = 8.The scattered data is measured at the point (x N , ζ), where ζ is a fixed height.The grid points for generating scattered data and surface reconstruction are set different.The reconstruction domain is discretized to N subintervals at all time steps.Noisy data is also used to test the robustness of the method: If the scattered field measured by the receiver is ψ s , then the noisy data is given by where ϵ is the noise level and ϑ is the random number in [−1, 1] generated with Gaussian distribution.
The rough surfaces used in the tests are randomly generated by a given autocorrelation function (a.c.f.) ρ(η), where η = x − x ′ .Two types of rough surfaces are employed here, having Gaussian-type a.c.f.where σ is the variance and l is the surface scale (autocorrelation length); and 'sub-fractal' a.c.f.
For both cases, we set σ = 0.2.The typical peak-to-trough of the tested rough surfaces is around 0.6, around λ/10.Qualitatively, at fine scales, the Gaussian surface is relatively smooth, whereas the sub-fractal surface is more jagged.Within the reconstruction algorithm, all functions are initially allowed to be complexvalued.As expected it is found that the reconstructed surface has only a small imaginary component, which can be neglected.The discretization size N is varied to test the performance of the method.The surface moves in the negative x-direction by ∆x each step for ∆x = L/N .The ℓ 2 -norm is used to quantify and examine the error, with where H is the reconstructed surface and Ĥ is the exact surface.Note that this can underestimate the qualitative performance, for example it may exaggerate the error due to a small horizontal shift in the reconstructed surface.

Dirichlet boundary condition
We first assume Dirichlet boundary conditions.The discretization size is chosen to be N = 500.The receiver is placed at height ζ = 0.7, which is about 0.11λ.The horizontal surface scale of the roughness l is initially set as the same as the wavelength, although this will later be varied.Initially the measurements are assumed to be noise-free.In the absence of noise the reconstructed surface against the exact surface for both types of surfaces is plotted in fig. 3. Clearly, the reconstructed surface agrees closely with the exact surface.A small height discrepancy can be seen between the reconstruction and exact surface near the peaks, but all key features are well-captured and correctly located horizontally.Table 1 gives the ℓ 2 -norm error between the actual and recovered surfaces with different discretization sizes N and different receiver heights ζ for the fixed surface scale l = λ.Clearly, with more nodes used, the algorithm is more effective.The receiver at height Gaussian type surface sub-fractal type surface reconstruction domain, which leads to a mild increase in error.
In order to test the robustness of the approach, we now add 3% noise to the scattered data at each time step, and the resulting amplitude of scattered field with and without noise is shown in fig. 4 using N = 500, ζ = 0.7, and l = λ for a Gaussian type surface.With the noisy data, the recovered surfaces compared to the actual surfaces are shown in fig. 5.With the noisy data, the reconstructed surfaces capture all main features, but before further processing display fine scale fluctuations throughout the domain, which for the most part are qualitatively similar to the measurement noise.However a few much larger spikes appear.Their presence is due to the ambiguity of taking the logarithm in eq. ( 19).On the other hand, it is observed that there are only a small number of these spikes in each case.To eliminate this problem, we introduce the following post-reconstruction 'filter' into the algorithm: Suppose that the surface heights recovered at the time steps t j−1 and t j are H j−1 and H j obtained by eqs.(19) and (20).We force the surface height H j equal to H j−1 if the surface derivative at this point is larger than one, namely |(H j − H j−1 )/∆x| > 1.With this filter, the reconstructed surfaces are shown in fig.6.A straightforward FFTtype filter can be applied to damp out most oscillations.After filtering, the reconstructed surface recaptures the overall shape of the original surface.We also list the ℓ 2 -norm error between the filtered recovered surfaces and the actual surfaces in table 3. The ℓ 2 -norm error increases as the noise level goes up.The error is larger for the recovered sub-fractal surface.

Neumann boundary condition
The applicability of the approach for Neumann boundary condition is now examined.With the discretization size N = 500, surface scale l = λ, and the receiver height ζ = 0.7, the     reconstructed surface versus the true surface is plotted in fig.7 for Gaussian and sub-fractal type surfaces under the Neumann boundary condition, in the absence of measurement noise.Similar to the Dirichlet case, the recovered surface closely matches the original surface, with the largest errors apparent around the peaks.
Fix the surface scale l = λ, the ℓ 2 -norm error for using different sizes of discretization and heights of receiver for both types of surfaces is summarised in table 4. The size of Gaussian type surface sub-fractal type surface  A white noise is added to the scattered data at each time step with the noise level 3%, the magnitude of the data with and without noise is shown in fig.8 with N = 500, ζ = 0.7 and l = λ for sub-fractal type surface.The reconstruction results for both type surfaces with noisy data are shown in fig.9.The method for the Neumann boundary condition is more robust with respect to the noise compared to the Dirichlet case.Without a filter, the reconstructed surface, though contains certain oscillations, still closely matches the actual surface.For the sub-fractal case, while the solution clearly captures the main surface peaks and troughs, the noise initially masks the smaller scale features.However, with a straightforward FFT type filter, the reconstruction can capture most shapes of the true surface.Finally, the ℓ 2 -norm error between the filtered rough surface and the actual surface of two different types with the noise data is listed in table 6  increases more, the increment in the error is clear.

Conclusions
A novel reconstruction algorithm has been developed to recover the rough surface based on repeated measurement of wave field at a single receiver.The surface is moved with respect to source and receiver during the algorithm, and the method is of a marchingtype, in which the surface height is recovered progressively as the position changes.This approach is based on the parabolic wave integral equations.Both Dirichlet and Neumann boundary conditions have been tackled.The solution relies on low grazing angle.Although the parabolic equation Green's function which neglects backscattering, a similar approach using the full Green's function is feasible.
Through numerical examples, it was found that the solution surface captures most shapes of the original surface for both Gaussian and sub-fractal type surfaces.Good agreement between recovered and real surface was observed, whereas large error shows up around peaks of surface.The results with respect to random white noise was also tested.Though some instabilities are present in the recovered surface under the Dirichlet boundary condition, the large peak-type errors can be damped out completely by employing a derivative filter in the algorithm.The noise induced in the solution is qualitatively similar to the measurement noise, i.e. it is delta-correlated.With a FFT-type filter, the reconstructed surface closely matches the exact solution.
A main motivation for this work is that the principles of the method can be extended in two important directions.Work is currently underway to adapt the approach to phaseless data.If the measurement receiver gives the amplitude of the wave field the reconstruction algorithm may be far more challenging.Some related previous work can be found in [30], but in which the surface is stationary and large number of receivers are required.In addition, the principles of this method are feasibly applicable to three-dimensional (3D) problems.In 3D, applying the standard marching method over a square N × N domain (say) would require O(N 2 ) receivers to obtain sufficient scattered data.Such a receiver array would in general be both difficult to accomplish and impractical without selfinterference.The approach here would replace this with a source and line array, moving with respect to the surface.At each step of the marching iteration the surface would be recovered along a line, orthogonal to the direction of the Gaussian beam.

Figure 2 :
Figure2: Illustration of reconstruction domain at successive time steps lying between source (S) and receiver (R).At time step t j , the surface in the reconstruction domain is h j (x).
ζ N ) at the fixed point r ζ N = (x N , ζ) where ζ is the height of the receiver for time steps t j , j = 1, 2, • • • , N .The algorithm seeks to reconstruct the unknown values to the right successively as the surface moves.

Figure 3 :
Figure 3: Reconstruction of rough surfaces with Dirichlet boundary condition for surfaces with (a) Gaussian type and (b) sub-fractal autocorrelation functions.

Figure 4 :Figure 5 :
Figure 4: Amplitude of scattered data at each time step without noise (upper) and with 3% noise (lower) under Dirichlet boundary condition. ℓ

Figure 6 :Figure 7 :
Figure 6: Reconstruction of rough surfaces under Dirichlet boundary condition with 3% scattered data using the derivative filter in the algorithm for (a) Gaussian type surface and (b) sub-fractal type surface; plus an additional FFT-type filter for (c) Gaussian type surface and (d) sub-fractal type surface.

Figure 8 :
Figure 8: Amplitude of scattered data at each time step without noise (upper) and with 3% noise (lower) under Neumann boundary condition.

Figure 9 :
Figure 9: Reconstruction of rough surfaces under Neumann boundary condition with 3% noisy scattered data for (a) Gaussian type surface and (b) sub-fractal type surface; and filtered recovered surfaces of (c) Gaussian type and (d) sub-fractal type using a FFT type filter.

Table 1 :
The ℓ 2 norm error between the actual and recovered surfaces for the Dirichlet boundary condition with different number of discretized points N and different heights of the receiver ζ. ζ = 0.7 gives the best result.The error is slightly larger for the sub-fractal type surface.The ℓ 2 -norm error for different values of surface scales (l) is listed in table 2 with fixed N = 500 and ζ = 0.7.As the surface scale decreases, more peaks occur within the

Table 2 :
The ℓ 2 -norm error between the actual surface and the recovered surface obtained using different surface scale sizes l and fixing N = 500 and ζ = 0.7 under the Dirichlet boundary condition.

Table 3 :
The ℓ 2 -norm error between the actual surface and the filtered recovered surface obtained by suing N = 500, ζ = 0.7 and l = λ under the Dirichlet boundary condition with noisy scattered data of different noise levels.

Table 4 :
The ℓ 2 norm error between the actual and recovered surfaces under the Neumann boundary condition with different number of discretized points N and different height of the receiver ζ.error is similar to the Dirichlet case.It is found that the method is robust with respect to the discretization size.Within this range, and with a larger number of points used, the error reduces.In the Neumann case, the error for both recovered Gaussian type and sub-fractal surfaces keeps at the similar level.The heights of the receiver at ζ = 0.7 and ζ = 1.0 give similar results while the error increases for ζ = 1.2.Now, we set N = 500 and ζ = 0.7, the error with respect to different values of surface scales (l) is listed in table 5. Unsurprisingly, the performance degrades as the surface

Table 5 :
The ℓ 2 -norm error between the actual and the recovered surfaces obtained by using different surface scale sizes l and fixing N = 500 and ζ = 0.7 under the Neumann boundary condition.

Table 6 :
. If the noise level The ℓ 2 -norm error between the actual surface and the filtered recovered surface obtained by suing N = 500,ζ = 0.7 and l = λ under the Neumann boundary condition with noisy data of different noise levels.
is not big enough (≥ 5%), the error of recovered surface using noisy data keeps at the same level with non-noisy data under the Neumann boundary condition.As noise level