- Email: [email protected]
Reconstruction of a fractal rough surface
Physica D 213 (2006) 25–30 www.elsevier.com/locate/physd
Reconstruction of a fractal rough surface Zhijie Cai ∗ , Deqiang Chen, Shuai Lu School of Mathematical Sciences, Fudan University, Shanghai 200433, China Received 18 May 2005; received in revised form 12 October 2005; accepted 27 October 2005 Available online 28 November 2005 Communicated by S. Kai
Abstract Scattering from a rough surface has been widely used. In recent years, fractal geometry has been successfully applied to investigate the interaction between a fractal rough surface and a wave. In this paper, we discuss the reconstruction of a fractal rough surface with periodic boundary conditions. To reconstruct the fractal rough surface accurately, the fractal dimension is inverted directly. At the same time, an efficient direct searching method is given to avoid the difficulty of calculating the derivative of the minimal target function. The numerical simulations show that this method has high convergence speed and high precision. c 2005 Elsevier B.V. All rights reserved.
Keywords: Scattering; Fractal rough surface; Reconstruction; Inverse problem
1. Introduction During recent decades, scattering from a rough surface has been widely used, such as in the diagnosis of optical interfaces, radar cross section from sea and land surfaces, and acoustic return from underwater bottom. There have been several approaches to solving this problem. Early research of reconstructing grating profiles was based on measuring diffraction efficiency (see [1,2]). Wombell and DeSanto presented an inversion of a small-lope rough surface in the Kirchhoff approximation by using a plane wave at fixed frequency and incident angle (see [3,4]). Spivack used the parabolic approximated equation for Gaussian beam incidence at grazing angle and derived an inverse approach (see [5,6]). The surface profile of roughness was usually described by periodic or random functions, such as the sine function and Gaussian function. However, it is found that the elaborate structure characteristic of a radar echo has been destroyed by the classic mathematics–physics method, and little information can be extracted from the radar echo. It is also clear that the typical fractal characteristic always exists in the
elaborate structure. In recent years, fractal geometry has been successfully applied to investigate the interaction between a fractal rough surface and a wave, since it can describe a very complex object by using several fractal characteristic parameters. The fractal geometry serves as a bridge between a periodic function and a random function, and its geometrical profile is governed by several fractal parameters including the fractal dimension. Since the fractal geometry combines both long-range order and short-range disorder as the character of the rough surfaces, it is close to the real surfaces. In this paper, the minimal target function method is used to invert the fractal dimension which can be used to reconstruct the rough surface. Numerical results show that this method has high precision. This paper is organized as follows. In Section 2, the mathematical model of the problem is presented. In Section 3, the numerical algorithm for solving nonlinear optimal problems is given. In Section 4, the fractal dimension of the rough surface is inverted and the error analysis is given. The influence of some parameters on the inversion results is also considered. In Section 5, some conclusions of this paper are given. 2. Mathematical model of the problem
∗ Corresponding author. Tel.: +86 21 65642469.
E-mail address: [email protected] (Z. Cai). c 2005 Elsevier B.V. All rights reserved. 0167-2789/$ - see front matter doi:10.1016/j.physd.2005.10.011
In recent years, some researchers have pointed out that for real fractal objects, the measured fractal dimension will
26
Z. Cai et al. / Physica D 213 (2006) 25–30
modes and evanescent fields. Since the evanescent wave decays rapidly, only a few evanescent modes are needed. Suppose that the scattering wave from a rough surface creates 2n discrete Floquet modes and every mode has scattering angle θm = sin−1 Lλ m , m = ±1, ±2, . . . , ±n, then using the Euler formula, (2.2) can be written as 2π m 0 1 Ψ (x, z) = cos x − ki z (z − d) L 2πm 0 + i sin x − ki z (z − d) + a0 eikz0 (z−d) L n X 2π m 2π m x + bm sin x + am cos L L m=1 × eikzm (z−d) , Fig. 1. Geometrical structure of the problem. Region I is a homogeneous space; region II consists of discrete Floquet modes of scattering fields produced from the random and periodic rough surface.
(2.3)
where a0 , am , bm are the coefficients to be determined. In region II, the sum of some functions Ψ c , Ψ s is used to express the total field. Here Ψ c and Ψ s satisfy the following Helmholtz equations: (∇ 2 + k 2 )Ψmc (x, z) = 0,
m = 0, 1, 2, . . . , n,
(2.4)
m = 1, 2, . . . , n
(2.5)
vary with different yardsticks. In other words, if a different yardstick is adopted in the measurement, a different fractal dimension will be obtained (see [7]). Since the real fractal object does not have the property of infinite self-similar structure, before choosing the scale and fixing the critical point, particular analysis of the fractal object’s structure character (e.g. hierarchy) should be done. To avoid choosing a yardstick, the minimal target function method is used to invert the fractal dimension directly. Consider a plane wave incident on a one-dimensional fractal surface, which is described by the fractal function z = f (x) (x ∈ [− L2 , L2 ]) (see Fig. 1). z = z 0 is the position where we survey. The region can be divided into two parts: region I (z ≥ d, |x| ≤ L2 ) is a homogeneous space, and region II ( f (x) ≤ z ≤ d, |x| ≤ L2 ) consists of discrete Floquet modes of scattering fields produced from the rough surface. Consider a TE-polarized plane wave with an incident angle θi on a one-dimensional fractal rough conducting surface z = f (x). The wave equations are written as
where the superscripts c, s indicate the fields with cos 2πm L x and sin 2πm x , respectively. Then, the total field in region II L can be written as
ψi (x, z) = eiki x x−iki z (z−d) ,
Ψ 2 (x, z) = A0 Ψ0c (x, z) +
(2.1)
where ki x = k sin θi , ki z = k cos θi , k = 2π λ is the wave number, λ is the wavelength. In region I, the total field [8,9] is usually expressed as Ψ 1 (x, z) = ei
2πm 0 L
+
x−iki z (z−d)
n X
Bm−m 0 ei
2πm L x+ik zm (z−d)
,
(2.2)
m=−n
where the first term on the right-hand side is the incident field, and it can be chosen m 0 = Lλ cos θi ; the second term is the scattered field, Bm−m 0 are the coefficients to be determined. The truncation number n is determined by the propagation
(∇ + k 2
2
)Ψms (x, z)
= 0,
with the Dirichlet and periodic boundary conditions Ψmc(s) (x, f (x)) = 0, L c(s) c(s) L ,z , − , z = Ψm Ψm 2 2 ∂ c(s) L ∂ c(s) L Ψm − ,z = Ψm ,z , ∂x 2 ∂x 2 2π m Ψmc (x, d) = cos x , L 2πm Ψms (x, d) = sin x , L
n X
(2.6) (2.7) (2.8) (2.9) (2.10)
[Am Ψmc (x, z) + Bm Ψms (x, z)].
m=1
(2.11) Thus, using the FEM, the Helmholtz equations (2.4)–(2.10) can be solved to get Ψ 2 . At z = d, the interface between regions I and II, Ψ 1 and 2 Ψ satisfy the continuous boundary conditions Ψ 1 (x, d) = ∂ ∂ Ψ 2 (x, d), ∂z Ψ 1 (x, d) = ∂z Ψ 2 (x, d). From the continuous boundary conditions, we can obtain δmm 0 + am = Am ,
m = 0, 1, . . . , n,
(2.12)
iδmm 0 + bm = Bm ,
m = 1, . . . , n
(2.13)
and
Z. Cai et al. / Physica D 213 (2006) 25–30
Z ∂Ψ c 2 − δm0 L/2 A0 0 (x, d) −iki z δmm 0 + ik zm am = L ∂z −L/2 2πm 2 − δm0 × cos x dx + L L Z n c s L/2 X ∂Ψ M ∂Ψ M × AM (x, d) + B M (x, d) ∂z ∂z M=1 −L/2 2π m x dx, m = 0, 1, . . . , n, × cos L Z ∂Ψ c 2 L/2 ki z δmm 0 + ik zm bm = A0 0 (x, d) ∂z L −L/2 2π m 2 x dx + × sin L L Z n c s L/2 X ∂Ψ M ∂Ψ M × AM (x, d) + B M (x, d) ∂z ∂z −L/2 M=1 2πm x dx, m = 1, . . . , n, × sin L
−i
2 L
Z
27
L/2
∂Ψms 0
−L/2
∂z
(x, d) sin
2π m x (1 − δ0m 0 )dx, L
m = 1, . . . , n.
(2.14)
(2.17)
Thus, the coefficients am and bm can be determined and the total field Ψ 1 can be obtained (see [9]). To reconstruct the fractal rough surface, the minimal target function method is used to invert the fractal dimension. The minimal target function is usually written as min E(D) =
N X
(Ii (D) − Ii∗ )2 ,
(2.18)
i=1
(2.15)
where δmm 0 is the Kronecker delta. Substituting (2.12) and (2.13) into (2.14) and (2.15), the following linear equations of am and bm can be obtained:
where E(D) is the square error function, D is the fractal dimension to be inverted, Ii (D) denote the intensities of the total field at the survey position z = z 0 , which can be calculated from Ψ 1 , Ii∗ are the measurements at z = z 0 , and N are the number of measurements. Note that the function f (x) in the boundary condition is related to the fractal dimension D (cf. (4.1) in Section 4), so the intensities Ii are generally the functions of D. 3. The nonlinear optimal method
∂Ψ0c
2π m 2 − δm0 (x, d) cos x dx −ik zm am + a0 L ∂z L −L/2 " Z n c X 2 − δm0 L/2 ∂Ψ M 2π m aM + (x, d) cos x dx L L −L/2 ∂z M=1 # Z s 2π m 2 − δm0 L/2 ∂Ψ M (x, d) cos x dx + bM L L −L/2 ∂z Z 2 − δm0 L/2 ∂Ψ0c (x, d) = −iki z δmm 0 − δ0m 0 L −L/2 ∂z Z 2πm 2 − δm0 L/2 ∂Ψmc 0 × cos x dx − (x, d) L L −L/2 ∂z 2π m 0 × cos x (1 − δ0m 0 )dx L Z 2 − δm0 L/2 ∂Ψms 0 −i (x, d) L −L/2 ∂z 2π m 0 × cos x (1 − δ0m 0 )dx, m = 0, 1, . . . , n, (2.16) L Z 2 L/2 ∂Ψ0c 2πm − ik zm bm + a0 (x, d) sin x dx L −L/2 ∂z L " Z n c X 2 L/2 ∂Ψ M 2π m + aM (x, d) sin x dx L −L/2 ∂z L M=1 # Z L/2 s ∂Ψ M 2 2πm + bM (x, d) sin x dx L −L/2 ∂z L Z 2 L/2 ∂Ψ0c 2πm = ki z δmm 0 − δ0m 0 (x, d) sin x dx L −L/2 ∂z L Z 2 L/2 ∂Ψmc 0 2π m 0 − (x, d) sin x (1 − δ0m 0 )dx L −L/2 ∂z L Z
L/2
In Section 2, the inverse problem of the fractal dimension is transformed into the optimal problem of the square error function. Generally, Ii (D) are calculated from a complicated partial differential equation. It is very difficult to calculate Ii0 (D) and large errors are usually caused by the approximate algorithms of the derivative. Therefore, it is necessary to adopt some direct searching methods that do not contain the differential information. Consider the following general optimal problem. Let f (x) be a continuous function defined in [a, b], to find the minimum point xmin , such that f (xmin ) = min f (x). x∈[a,b]
(3.1)
The direct searching algorithm is as follows. Step 1. Choose the initial data x10 , x20 ∈ [a, b], x10 6= x20 . Without loss of generality, assume that f (x10 ) ≤ f (x20 ). Let the iterative step k = 0. Step 2. If |x1k − x2k | < ε (ε is a given precision), go to Step 8. Step 3. Let x∗k = x1k − η(x2k − x1k ). Step 4. If f (x∗k ) ≤ f (x1k ) and x∗k ∈ [a, b], then let x1k+1 = x∗k , x2k+1 = x1k , k = k + 1, go to Step 2. k Step 5. If f (x∗k ) > f (x1k ) or x∗k 6∈ [a, b], let x∗∗ = k k k x1 + η(x2 − x1 ). k ) < f (x k ), then let x k+1 = x k , x k+1 = x k , Step 6. If f (x∗∗ ∗∗ 2 1 1 1 k = k + 1, go to Step 2. k ) < f (x k ), then let x k+1 = x k , Step 7. If f (x1k ) ≤ f (x∗∗ 2 1 1 k , k = k + 1, go to Step 2. x2k+1 = x∗∗ Step 8. Let x ∗ = x1k , stop the algorithm. The following theorem ensures the convergence of the algorithm.
28
Z. Cai et al. / Physica D 213 (2006) 25–30
(a) D = 1.3.
(b) D = 1.5. Fig. 2. The curves of the fractal function with different fractal dimensions.
Theorem 3.1 ([10]). Assume that f (x) is a continuous function defined in [a, b]. It has a unique minimum point xmin , and the corresponding minimal value is f min = f (xmin ). If the parameter η satisfies ! |x20 − x10 | 1 , ≤ η < 1, (3.2) max 1 − b−a 2 then the sequences {x1k }, {x2k } generated by the above algorithm converge to the minimum point xmin . To find the global minimum of the inverse problem (2.18), the Monte Carlo method is used to generate several pairs of initial data D10 , D20 to solve the inverse problem. Then the minimal one is chosen as the approximate value of the global minimum. 4. Numerical results 4.1. Minimal target function There are several functions describing the fractal surface. The band-limited Weierstrass–Mandelbrot function [11–14] is chosen in this paper: f (x) = δC
Real dimension Initial data Simulation result Iterative step RMSE Maximum error
D = 1.3 D10 = 1.25, D20 = 1.26 D ∗ = 1.2998 26 1.0083 × 10−5 2.3988 × 10−5
D = 1.5 D10 = 1.53, D20 = 1.535 D ∗ = 1.4997 31 1.7522 × 10−5 5.3062 × 10−5
This means that the curve looks similar to the original one when the horizontal axis is scaled by b and the vertical axis 1 by D−1 . Fig. 2 shows the curves of the fractal function with different fractal dimensions. With the fractal dimension increasing, it is found that the curve becomes more and more complex. For a real fractal object, the measured fractal dimension will vary with different yardsticks. To avoid choosing a yardstick, the minimal target function method is adopted to invert the fractal dimension directly. As mentioned in Section 2, the least square error function (2.18) is chosen as the minimal target function. 4.2. Simulation results of the fractal dimension
M−1 X
(D − 1) sin(k0 b x + ψm ), m
m
(4.1)
m=0
where D (1 < D < 2) is the fractal dimension, ψm is the random phase, b (b > 1) is the spatial frequency scaling parameter, k0 is the fundamental spatial wavep number, and δ is the rms surface height. The coefficient C = 2D(2 − D)/(1 − (D − 1)2M ) is a normalized amplifier control parameter, and M is the highest harmonic wave number. The fractal function has a finite band of spatial frequency, and owns self-similarity in the corresponding finite range of differentiate rates. With increasing the frequency of the periodic function, the elaborate structure is described using a weighted sum of periodic functions. If ψm = 0, the self-similarity of (4.1) is reflected by the following function: f (bx) ≈
Table 1 Comparison between the calculated dimensions and the real measurements
1 f (x). D−1
(4.2)
Using the above algorithm, some simulations are made with D = 1.3 and D = 1.5, respectively. The comparison between the inverted dimension and the real one of the fractal rough surface is shown in Table 1. The root mean square error (RMSE) and maximum error between the reconstructed surface and the original real one are also shown in Table 1. Table 1 shows that the inverted dimension is close to the real one. Meanwhile, the intensity of the total field and the rough surface calculated from the inverted dimension also accord with that calculated from the corresponding real one. The errors of the rough surface are shown in Figs. 3 and 4. From these two figures we can see that the errors between the simulated and the real fractal rough surface are very small. It is found that our approach can reconstruct the rough surface very well.
Z. Cai et al. / Physica D 213 (2006) 25–30
Fig. 3. The error of the fractal rough surface between D = 1.3 and D ∗ = 1.2998.
Fig. 4. The error of the fractal rough surface between D = 1.5 and D ∗ = 1.4997. Table 2 Simulation results with different η Real dimension Result of η = 0.5 RMSE Maximum error Result of η = 0.8 RMSE Maximum error
D = 1.3 1.2998 1.0083 × 10−5 2.3988 × 10−5 1.2998 1.0083 × 10−5 2.3988 × 10−5
D = 1.5 1.4997 1.7522 × 10−5 5.3062 × 10−5 1.5003 1.7528 × 10−5 5.3092 × 10−5
4.3. The influence of the parameters In this subsection, we analyze the influence of some parameters, such as η, n, z 0 , and the error of the measurements. 4.3.1. The influence of η on the simulation results Select η = 0.5 and η = 0.8 in the algorithm, respectively. The inversion results are shown in Table 2. From Table 2, it is found that the inversion results with different η are very close. It is easy to prove that the
29
Fig. 5. The error of the fractal rough surface between D = 1.5 and D ∗ = 1.5003.
Fig. 6. The error of the fractal rough surface between D = 1.3 and D ∗ = 1.3001.
convergence rate of the algorithm is − logη b−a ε . So the smaller the η we choose, the faster the convergence rate of the method. From Theorem 3.1, η ∈ [ 12 , 1) can be chosen arbitrarily, so D10 = 1, D20 = 2, η = 0.5 are chosen in the following calculations. Figs. 4 and 5 show the error of the rough surface between D = 1.5 and D ∗ = 1.4997, 1.5003, respectively. 4.3.2. The influence of n on the simulation results Here n is the number of expanding coefficients of the total field. To study the influence of n, take n = 30 and n = 20, respectively. The calculation results are shown in Table 3. From Table 3, it is found that the number of expanding coefficients of the total field n has little influence on the simulation results. Therefore, a small value of n can be used in the calculation to reduce the computing time. The errors of the rough surface between D = 1.3 and D ∗ = 1.2998, 1.3001 are shown in Figs. 3 and 6.
30
Z. Cai et al. / Physica D 213 (2006) 25–30
In the 2-D case, 2-D Weierstrass–Mandelbrot functions are usually used to describe the 2-D fractal surface:
Table 3 Simulation results with different n Real dimension Result of n = 30 RMSE Maximum error Result of n = 20 RMSE Maximum error
D = 1.3 1.2998 1.0083 × 10−5 2.3988 × 10−5 1.3001 5.0420 × 10−6 1.1997 × 10−5
D = 1.5 1.4997 1.7522 × 10−5 5.3062 × 10−5 1.5003 1.7528 × 10−5 5.3092 × 10−5
Table 4 Simulation results with different survey positions z 0 Real dimension Result of z 0 = 18λ RMSE Maximum error Result of z 0 = 1000λ RMSE Maximum error
D = 1.3 1.2998 1.0083 × 10−5 2.3988 × 10−5 1.3003 1.5127 × 10−5 3.5999 × 10−5
D = 1.5 1.4997 1.7522 × 10−5 5.3062 × 10−5 1.5001 5.8421 × 10−6 1.7694 × 10−5
Table 5 The influence of measurements error on the simulation results (D = 1.3) Error level (%)
Simulation result
RMSE
Maximum error
0.5 1 2 3 5
1.2998 1.3004 1.2995 1.3183 1.3162
1.0083 × 10−5 2.0169 × 10−5 2.5207 × 10−5 9.2631 × 10−4 8.1963 × 10−4
2.3988 × 10−5 4.8004 × 10−5 5.9951 × 10−5 2.2355 × 10−3 1.9748 × 10−3
4.3.3. The influence of survey position z 0 on the simulation results In studying the influence of survey position z 0 , let z 0 = 18λ and z 0 = 1000λ, respectively. The inversion results are shown in Table 4. It is found that the survey position also has little influence on the simulation results. 4.3.4. The influence of measurements error on the simulation results To test the influence of measurements error, some random errors are added on the measurements. Table 5 shows the results with different random error levels. From Table 5, the maximum error is less than 1% when the measurements error is 5%; the results are quite satisfactory. 5. Conclusions In this paper, the minimal target function method is used to invert the fractal dimension of a fractal rough surface directly, and then the rough surface is reconstructed. Numerical simulations show that this method can achieve high accuracy, and the optimal algorithm is very effective; the derivative of the target function need not be calculated. Also, the measurements error and the parameters such as η, n and z 0 have little influence on the simulation results. So certain parameters can be chosen freely in practice.
f (x, y) = δC
M−1 X
(D − 1)m cos[k0 bm (x + y) + ψm ],
(5.1)
m=0
where D (1 < D < 2) is the fractal dimension, ψm is the random phase, b (b > 1) is the spatial frequency scaling parameter, k0 is the fundamental spatial wavep number, and δ is the rms surface height. The coefficient C = 2D(2 − D)/(1 − (D − 1)2M ) is a normalized amplifier control parameter, and M is the highest harmonic wave number. Supposing all the parameters except D are known, the fractal dimension D can be inverted using our approach, as we did in the case of the 1-D surface profile. Acknowledgement The project is supported by the National Natural Science Foundation of China (10201008). References [1] A. Roger, D. Maystre, Inverse scattering method in electromagnetic optics: application to diffraction gratings, J. Opt. Soc. Amer. 70 (12) (1980) 1483–1495. [2] A.M. Husier, A. Quatropani, H.D. Baltes, Construction of grating profiles yielding prescribed diffraction efficiency, Opt. Commun. 41 (1982) 149–153. [3] R.J. Wombell, J.A. DeSanto, The reconstruction of shallow rough-surface profiles from scattered field data, Inverse Problems 7 (1) (1991) L7–L12. [4] R.J. Wombell, J.A. DeSanto, Reconstruction of rough-surface profiles with the Kirchhoff approximation, J. Opt. Soc. Amer. A 8 (12) (1991) 1892–1897. [5] M. Spivack, Direct solution of the inverse problem for rough surface scattering at grazing incidence, J. Phys. A 25 (11) (1992) 3295–3302. [6] M. Spivack, A numerical approach to rough-surface scattering by parabolic equation method, J. Acoust. Soc. Amer. 87 (5) (1990) 1999–2004. [7] M. Spivack, Solution of the inverse-scattering problem for grazing incidence upon a rough surface, J. Opt. Soc. Amer. A 9 (8) (1992) 1352–1355. [8] Y.Q. Jin, G. Li, Detection of a scatter target over a randomly rough surface by using the angular correlation function in a finite-element approach, Waves Random Media 10 (2000) 273–280. [9] S.H. Luo, L. Tsang, C.H. Chan, Application of the finite element method to Monte Carlo simulations of scattering of waves by random rough surfaces with the periodic condition, J. Electromagn. Waves Appl. 5 (8) (1991) 835–855. [10] Z. Cai, D. Chen, A direct searching method for nonlinear optimal problems and its convergence proof (in Chinese) (submitted for publication). [11] D.L. Jaggard, X. Sun, Scattering from fractally corrugated surfaces, J. Opt. Soc. Amer. A 7 (6) (1990) 1131–1139. [12] Z. Li, Y. Jin, Numerical simulation of bistatic scattering from fractal rough surface in the finite element method, Sci. China Ser. E 44 (1) (2001) 12–18. [13] Y.Q. Jin, Z. Li, Bistatic scattering and transmitting through fractal rough dielectric surface using FBM/SAA method, J. Electromagn. Waves Appl. 16 (4) (2002) 551–572. [14] Y.Q. Jin, Z. Li, Reconstruction of roughness profile of fractal surfaces from scattering measurement at grazing incidence, J. Appl. Phys. 89 (3) (2001) 1922–1926.
Reconstruction of a fractal rough surface Zhijie Cai ∗ , Deqiang Chen, Shuai Lu School of Mathematical Sciences, Fudan University, Shanghai 200433, China Received 18 May 2005; received in revised form 12 October 2005; accepted 27 October 2005 Available online 28 November 2005 Communicated by S. Kai
Abstract Scattering from a rough surface has been widely used. In recent years, fractal geometry has been successfully applied to investigate the interaction between a fractal rough surface and a wave. In this paper, we discuss the reconstruction of a fractal rough surface with periodic boundary conditions. To reconstruct the fractal rough surface accurately, the fractal dimension is inverted directly. At the same time, an efficient direct searching method is given to avoid the difficulty of calculating the derivative of the minimal target function. The numerical simulations show that this method has high convergence speed and high precision. c 2005 Elsevier B.V. All rights reserved.
Keywords: Scattering; Fractal rough surface; Reconstruction; Inverse problem
1. Introduction During recent decades, scattering from a rough surface has been widely used, such as in the diagnosis of optical interfaces, radar cross section from sea and land surfaces, and acoustic return from underwater bottom. There have been several approaches to solving this problem. Early research of reconstructing grating profiles was based on measuring diffraction efficiency (see [1,2]). Wombell and DeSanto presented an inversion of a small-lope rough surface in the Kirchhoff approximation by using a plane wave at fixed frequency and incident angle (see [3,4]). Spivack used the parabolic approximated equation for Gaussian beam incidence at grazing angle and derived an inverse approach (see [5,6]). The surface profile of roughness was usually described by periodic or random functions, such as the sine function and Gaussian function. However, it is found that the elaborate structure characteristic of a radar echo has been destroyed by the classic mathematics–physics method, and little information can be extracted from the radar echo. It is also clear that the typical fractal characteristic always exists in the
elaborate structure. In recent years, fractal geometry has been successfully applied to investigate the interaction between a fractal rough surface and a wave, since it can describe a very complex object by using several fractal characteristic parameters. The fractal geometry serves as a bridge between a periodic function and a random function, and its geometrical profile is governed by several fractal parameters including the fractal dimension. Since the fractal geometry combines both long-range order and short-range disorder as the character of the rough surfaces, it is close to the real surfaces. In this paper, the minimal target function method is used to invert the fractal dimension which can be used to reconstruct the rough surface. Numerical results show that this method has high precision. This paper is organized as follows. In Section 2, the mathematical model of the problem is presented. In Section 3, the numerical algorithm for solving nonlinear optimal problems is given. In Section 4, the fractal dimension of the rough surface is inverted and the error analysis is given. The influence of some parameters on the inversion results is also considered. In Section 5, some conclusions of this paper are given. 2. Mathematical model of the problem
∗ Corresponding author. Tel.: +86 21 65642469.
E-mail address: [email protected] (Z. Cai). c 2005 Elsevier B.V. All rights reserved. 0167-2789/$ - see front matter doi:10.1016/j.physd.2005.10.011
In recent years, some researchers have pointed out that for real fractal objects, the measured fractal dimension will
26
Z. Cai et al. / Physica D 213 (2006) 25–30
modes and evanescent fields. Since the evanescent wave decays rapidly, only a few evanescent modes are needed. Suppose that the scattering wave from a rough surface creates 2n discrete Floquet modes and every mode has scattering angle θm = sin−1 Lλ m , m = ±1, ±2, . . . , ±n, then using the Euler formula, (2.2) can be written as 2π m 0 1 Ψ (x, z) = cos x − ki z (z − d) L 2πm 0 + i sin x − ki z (z − d) + a0 eikz0 (z−d) L n X 2π m 2π m x + bm sin x + am cos L L m=1 × eikzm (z−d) , Fig. 1. Geometrical structure of the problem. Region I is a homogeneous space; region II consists of discrete Floquet modes of scattering fields produced from the random and periodic rough surface.
(2.3)
where a0 , am , bm are the coefficients to be determined. In region II, the sum of some functions Ψ c , Ψ s is used to express the total field. Here Ψ c and Ψ s satisfy the following Helmholtz equations: (∇ 2 + k 2 )Ψmc (x, z) = 0,
m = 0, 1, 2, . . . , n,
(2.4)
m = 1, 2, . . . , n
(2.5)
vary with different yardsticks. In other words, if a different yardstick is adopted in the measurement, a different fractal dimension will be obtained (see [7]). Since the real fractal object does not have the property of infinite self-similar structure, before choosing the scale and fixing the critical point, particular analysis of the fractal object’s structure character (e.g. hierarchy) should be done. To avoid choosing a yardstick, the minimal target function method is used to invert the fractal dimension directly. Consider a plane wave incident on a one-dimensional fractal surface, which is described by the fractal function z = f (x) (x ∈ [− L2 , L2 ]) (see Fig. 1). z = z 0 is the position where we survey. The region can be divided into two parts: region I (z ≥ d, |x| ≤ L2 ) is a homogeneous space, and region II ( f (x) ≤ z ≤ d, |x| ≤ L2 ) consists of discrete Floquet modes of scattering fields produced from the rough surface. Consider a TE-polarized plane wave with an incident angle θi on a one-dimensional fractal rough conducting surface z = f (x). The wave equations are written as
where the superscripts c, s indicate the fields with cos 2πm L x and sin 2πm x , respectively. Then, the total field in region II L can be written as
ψi (x, z) = eiki x x−iki z (z−d) ,
Ψ 2 (x, z) = A0 Ψ0c (x, z) +
(2.1)
where ki x = k sin θi , ki z = k cos θi , k = 2π λ is the wave number, λ is the wavelength. In region I, the total field [8,9] is usually expressed as Ψ 1 (x, z) = ei
2πm 0 L
+
x−iki z (z−d)
n X
Bm−m 0 ei
2πm L x+ik zm (z−d)
,
(2.2)
m=−n
where the first term on the right-hand side is the incident field, and it can be chosen m 0 = Lλ cos θi ; the second term is the scattered field, Bm−m 0 are the coefficients to be determined. The truncation number n is determined by the propagation
(∇ + k 2
2
)Ψms (x, z)
= 0,
with the Dirichlet and periodic boundary conditions Ψmc(s) (x, f (x)) = 0, L c(s) c(s) L ,z , − , z = Ψm Ψm 2 2 ∂ c(s) L ∂ c(s) L Ψm − ,z = Ψm ,z , ∂x 2 ∂x 2 2π m Ψmc (x, d) = cos x , L 2πm Ψms (x, d) = sin x , L
n X
(2.6) (2.7) (2.8) (2.9) (2.10)
[Am Ψmc (x, z) + Bm Ψms (x, z)].
m=1
(2.11) Thus, using the FEM, the Helmholtz equations (2.4)–(2.10) can be solved to get Ψ 2 . At z = d, the interface between regions I and II, Ψ 1 and 2 Ψ satisfy the continuous boundary conditions Ψ 1 (x, d) = ∂ ∂ Ψ 2 (x, d), ∂z Ψ 1 (x, d) = ∂z Ψ 2 (x, d). From the continuous boundary conditions, we can obtain δmm 0 + am = Am ,
m = 0, 1, . . . , n,
(2.12)
iδmm 0 + bm = Bm ,
m = 1, . . . , n
(2.13)
and
Z. Cai et al. / Physica D 213 (2006) 25–30
Z ∂Ψ c 2 − δm0 L/2 A0 0 (x, d) −iki z δmm 0 + ik zm am = L ∂z −L/2 2πm 2 − δm0 × cos x dx + L L Z n c s L/2 X ∂Ψ M ∂Ψ M × AM (x, d) + B M (x, d) ∂z ∂z M=1 −L/2 2π m x dx, m = 0, 1, . . . , n, × cos L Z ∂Ψ c 2 L/2 ki z δmm 0 + ik zm bm = A0 0 (x, d) ∂z L −L/2 2π m 2 x dx + × sin L L Z n c s L/2 X ∂Ψ M ∂Ψ M × AM (x, d) + B M (x, d) ∂z ∂z −L/2 M=1 2πm x dx, m = 1, . . . , n, × sin L
−i
2 L
Z
27
L/2
∂Ψms 0
−L/2
∂z
(x, d) sin
2π m x (1 − δ0m 0 )dx, L
m = 1, . . . , n.
(2.14)
(2.17)
Thus, the coefficients am and bm can be determined and the total field Ψ 1 can be obtained (see [9]). To reconstruct the fractal rough surface, the minimal target function method is used to invert the fractal dimension. The minimal target function is usually written as min E(D) =
N X
(Ii (D) − Ii∗ )2 ,
(2.18)
i=1
(2.15)
where δmm 0 is the Kronecker delta. Substituting (2.12) and (2.13) into (2.14) and (2.15), the following linear equations of am and bm can be obtained:
where E(D) is the square error function, D is the fractal dimension to be inverted, Ii (D) denote the intensities of the total field at the survey position z = z 0 , which can be calculated from Ψ 1 , Ii∗ are the measurements at z = z 0 , and N are the number of measurements. Note that the function f (x) in the boundary condition is related to the fractal dimension D (cf. (4.1) in Section 4), so the intensities Ii are generally the functions of D. 3. The nonlinear optimal method
∂Ψ0c
2π m 2 − δm0 (x, d) cos x dx −ik zm am + a0 L ∂z L −L/2 " Z n c X 2 − δm0 L/2 ∂Ψ M 2π m aM + (x, d) cos x dx L L −L/2 ∂z M=1 # Z s 2π m 2 − δm0 L/2 ∂Ψ M (x, d) cos x dx + bM L L −L/2 ∂z Z 2 − δm0 L/2 ∂Ψ0c (x, d) = −iki z δmm 0 − δ0m 0 L −L/2 ∂z Z 2πm 2 − δm0 L/2 ∂Ψmc 0 × cos x dx − (x, d) L L −L/2 ∂z 2π m 0 × cos x (1 − δ0m 0 )dx L Z 2 − δm0 L/2 ∂Ψms 0 −i (x, d) L −L/2 ∂z 2π m 0 × cos x (1 − δ0m 0 )dx, m = 0, 1, . . . , n, (2.16) L Z 2 L/2 ∂Ψ0c 2πm − ik zm bm + a0 (x, d) sin x dx L −L/2 ∂z L " Z n c X 2 L/2 ∂Ψ M 2π m + aM (x, d) sin x dx L −L/2 ∂z L M=1 # Z L/2 s ∂Ψ M 2 2πm + bM (x, d) sin x dx L −L/2 ∂z L Z 2 L/2 ∂Ψ0c 2πm = ki z δmm 0 − δ0m 0 (x, d) sin x dx L −L/2 ∂z L Z 2 L/2 ∂Ψmc 0 2π m 0 − (x, d) sin x (1 − δ0m 0 )dx L −L/2 ∂z L Z
L/2
In Section 2, the inverse problem of the fractal dimension is transformed into the optimal problem of the square error function. Generally, Ii (D) are calculated from a complicated partial differential equation. It is very difficult to calculate Ii0 (D) and large errors are usually caused by the approximate algorithms of the derivative. Therefore, it is necessary to adopt some direct searching methods that do not contain the differential information. Consider the following general optimal problem. Let f (x) be a continuous function defined in [a, b], to find the minimum point xmin , such that f (xmin ) = min f (x). x∈[a,b]
(3.1)
The direct searching algorithm is as follows. Step 1. Choose the initial data x10 , x20 ∈ [a, b], x10 6= x20 . Without loss of generality, assume that f (x10 ) ≤ f (x20 ). Let the iterative step k = 0. Step 2. If |x1k − x2k | < ε (ε is a given precision), go to Step 8. Step 3. Let x∗k = x1k − η(x2k − x1k ). Step 4. If f (x∗k ) ≤ f (x1k ) and x∗k ∈ [a, b], then let x1k+1 = x∗k , x2k+1 = x1k , k = k + 1, go to Step 2. k Step 5. If f (x∗k ) > f (x1k ) or x∗k 6∈ [a, b], let x∗∗ = k k k x1 + η(x2 − x1 ). k ) < f (x k ), then let x k+1 = x k , x k+1 = x k , Step 6. If f (x∗∗ ∗∗ 2 1 1 1 k = k + 1, go to Step 2. k ) < f (x k ), then let x k+1 = x k , Step 7. If f (x1k ) ≤ f (x∗∗ 2 1 1 k , k = k + 1, go to Step 2. x2k+1 = x∗∗ Step 8. Let x ∗ = x1k , stop the algorithm. The following theorem ensures the convergence of the algorithm.
28
Z. Cai et al. / Physica D 213 (2006) 25–30
(a) D = 1.3.
(b) D = 1.5. Fig. 2. The curves of the fractal function with different fractal dimensions.
Theorem 3.1 ([10]). Assume that f (x) is a continuous function defined in [a, b]. It has a unique minimum point xmin , and the corresponding minimal value is f min = f (xmin ). If the parameter η satisfies ! |x20 − x10 | 1 , ≤ η < 1, (3.2) max 1 − b−a 2 then the sequences {x1k }, {x2k } generated by the above algorithm converge to the minimum point xmin . To find the global minimum of the inverse problem (2.18), the Monte Carlo method is used to generate several pairs of initial data D10 , D20 to solve the inverse problem. Then the minimal one is chosen as the approximate value of the global minimum. 4. Numerical results 4.1. Minimal target function There are several functions describing the fractal surface. The band-limited Weierstrass–Mandelbrot function [11–14] is chosen in this paper: f (x) = δC
Real dimension Initial data Simulation result Iterative step RMSE Maximum error
D = 1.3 D10 = 1.25, D20 = 1.26 D ∗ = 1.2998 26 1.0083 × 10−5 2.3988 × 10−5
D = 1.5 D10 = 1.53, D20 = 1.535 D ∗ = 1.4997 31 1.7522 × 10−5 5.3062 × 10−5
This means that the curve looks similar to the original one when the horizontal axis is scaled by b and the vertical axis 1 by D−1 . Fig. 2 shows the curves of the fractal function with different fractal dimensions. With the fractal dimension increasing, it is found that the curve becomes more and more complex. For a real fractal object, the measured fractal dimension will vary with different yardsticks. To avoid choosing a yardstick, the minimal target function method is adopted to invert the fractal dimension directly. As mentioned in Section 2, the least square error function (2.18) is chosen as the minimal target function. 4.2. Simulation results of the fractal dimension
M−1 X
(D − 1) sin(k0 b x + ψm ), m
m
(4.1)
m=0
where D (1 < D < 2) is the fractal dimension, ψm is the random phase, b (b > 1) is the spatial frequency scaling parameter, k0 is the fundamental spatial wavep number, and δ is the rms surface height. The coefficient C = 2D(2 − D)/(1 − (D − 1)2M ) is a normalized amplifier control parameter, and M is the highest harmonic wave number. The fractal function has a finite band of spatial frequency, and owns self-similarity in the corresponding finite range of differentiate rates. With increasing the frequency of the periodic function, the elaborate structure is described using a weighted sum of periodic functions. If ψm = 0, the self-similarity of (4.1) is reflected by the following function: f (bx) ≈
Table 1 Comparison between the calculated dimensions and the real measurements
1 f (x). D−1
(4.2)
Using the above algorithm, some simulations are made with D = 1.3 and D = 1.5, respectively. The comparison between the inverted dimension and the real one of the fractal rough surface is shown in Table 1. The root mean square error (RMSE) and maximum error between the reconstructed surface and the original real one are also shown in Table 1. Table 1 shows that the inverted dimension is close to the real one. Meanwhile, the intensity of the total field and the rough surface calculated from the inverted dimension also accord with that calculated from the corresponding real one. The errors of the rough surface are shown in Figs. 3 and 4. From these two figures we can see that the errors between the simulated and the real fractal rough surface are very small. It is found that our approach can reconstruct the rough surface very well.
Z. Cai et al. / Physica D 213 (2006) 25–30
Fig. 3. The error of the fractal rough surface between D = 1.3 and D ∗ = 1.2998.
Fig. 4. The error of the fractal rough surface between D = 1.5 and D ∗ = 1.4997. Table 2 Simulation results with different η Real dimension Result of η = 0.5 RMSE Maximum error Result of η = 0.8 RMSE Maximum error
D = 1.3 1.2998 1.0083 × 10−5 2.3988 × 10−5 1.2998 1.0083 × 10−5 2.3988 × 10−5
D = 1.5 1.4997 1.7522 × 10−5 5.3062 × 10−5 1.5003 1.7528 × 10−5 5.3092 × 10−5
4.3. The influence of the parameters In this subsection, we analyze the influence of some parameters, such as η, n, z 0 , and the error of the measurements. 4.3.1. The influence of η on the simulation results Select η = 0.5 and η = 0.8 in the algorithm, respectively. The inversion results are shown in Table 2. From Table 2, it is found that the inversion results with different η are very close. It is easy to prove that the
29
Fig. 5. The error of the fractal rough surface between D = 1.5 and D ∗ = 1.5003.
Fig. 6. The error of the fractal rough surface between D = 1.3 and D ∗ = 1.3001.
convergence rate of the algorithm is − logη b−a ε . So the smaller the η we choose, the faster the convergence rate of the method. From Theorem 3.1, η ∈ [ 12 , 1) can be chosen arbitrarily, so D10 = 1, D20 = 2, η = 0.5 are chosen in the following calculations. Figs. 4 and 5 show the error of the rough surface between D = 1.5 and D ∗ = 1.4997, 1.5003, respectively. 4.3.2. The influence of n on the simulation results Here n is the number of expanding coefficients of the total field. To study the influence of n, take n = 30 and n = 20, respectively. The calculation results are shown in Table 3. From Table 3, it is found that the number of expanding coefficients of the total field n has little influence on the simulation results. Therefore, a small value of n can be used in the calculation to reduce the computing time. The errors of the rough surface between D = 1.3 and D ∗ = 1.2998, 1.3001 are shown in Figs. 3 and 6.
30
Z. Cai et al. / Physica D 213 (2006) 25–30
In the 2-D case, 2-D Weierstrass–Mandelbrot functions are usually used to describe the 2-D fractal surface:
Table 3 Simulation results with different n Real dimension Result of n = 30 RMSE Maximum error Result of n = 20 RMSE Maximum error
D = 1.3 1.2998 1.0083 × 10−5 2.3988 × 10−5 1.3001 5.0420 × 10−6 1.1997 × 10−5
D = 1.5 1.4997 1.7522 × 10−5 5.3062 × 10−5 1.5003 1.7528 × 10−5 5.3092 × 10−5
Table 4 Simulation results with different survey positions z 0 Real dimension Result of z 0 = 18λ RMSE Maximum error Result of z 0 = 1000λ RMSE Maximum error
D = 1.3 1.2998 1.0083 × 10−5 2.3988 × 10−5 1.3003 1.5127 × 10−5 3.5999 × 10−5
D = 1.5 1.4997 1.7522 × 10−5 5.3062 × 10−5 1.5001 5.8421 × 10−6 1.7694 × 10−5
Table 5 The influence of measurements error on the simulation results (D = 1.3) Error level (%)
Simulation result
RMSE
Maximum error
0.5 1 2 3 5
1.2998 1.3004 1.2995 1.3183 1.3162
1.0083 × 10−5 2.0169 × 10−5 2.5207 × 10−5 9.2631 × 10−4 8.1963 × 10−4
2.3988 × 10−5 4.8004 × 10−5 5.9951 × 10−5 2.2355 × 10−3 1.9748 × 10−3
4.3.3. The influence of survey position z 0 on the simulation results In studying the influence of survey position z 0 , let z 0 = 18λ and z 0 = 1000λ, respectively. The inversion results are shown in Table 4. It is found that the survey position also has little influence on the simulation results. 4.3.4. The influence of measurements error on the simulation results To test the influence of measurements error, some random errors are added on the measurements. Table 5 shows the results with different random error levels. From Table 5, the maximum error is less than 1% when the measurements error is 5%; the results are quite satisfactory. 5. Conclusions In this paper, the minimal target function method is used to invert the fractal dimension of a fractal rough surface directly, and then the rough surface is reconstructed. Numerical simulations show that this method can achieve high accuracy, and the optimal algorithm is very effective; the derivative of the target function need not be calculated. Also, the measurements error and the parameters such as η, n and z 0 have little influence on the simulation results. So certain parameters can be chosen freely in practice.
f (x, y) = δC
M−1 X
(D − 1)m cos[k0 bm (x + y) + ψm ],
(5.1)
m=0
where D (1 < D < 2) is the fractal dimension, ψm is the random phase, b (b > 1) is the spatial frequency scaling parameter, k0 is the fundamental spatial wavep number, and δ is the rms surface height. The coefficient C = 2D(2 − D)/(1 − (D − 1)2M ) is a normalized amplifier control parameter, and M is the highest harmonic wave number. Supposing all the parameters except D are known, the fractal dimension D can be inverted using our approach, as we did in the case of the 1-D surface profile. Acknowledgement The project is supported by the National Natural Science Foundation of China (10201008). References [1] A. Roger, D. Maystre, Inverse scattering method in electromagnetic optics: application to diffraction gratings, J. Opt. Soc. Amer. 70 (12) (1980) 1483–1495. [2] A.M. Husier, A. Quatropani, H.D. Baltes, Construction of grating profiles yielding prescribed diffraction efficiency, Opt. Commun. 41 (1982) 149–153. [3] R.J. Wombell, J.A. DeSanto, The reconstruction of shallow rough-surface profiles from scattered field data, Inverse Problems 7 (1) (1991) L7–L12. [4] R.J. Wombell, J.A. DeSanto, Reconstruction of rough-surface profiles with the Kirchhoff approximation, J. Opt. Soc. Amer. A 8 (12) (1991) 1892–1897. [5] M. Spivack, Direct solution of the inverse problem for rough surface scattering at grazing incidence, J. Phys. A 25 (11) (1992) 3295–3302. [6] M. Spivack, A numerical approach to rough-surface scattering by parabolic equation method, J. Acoust. Soc. Amer. 87 (5) (1990) 1999–2004. [7] M. Spivack, Solution of the inverse-scattering problem for grazing incidence upon a rough surface, J. Opt. Soc. Amer. A 9 (8) (1992) 1352–1355. [8] Y.Q. Jin, G. Li, Detection of a scatter target over a randomly rough surface by using the angular correlation function in a finite-element approach, Waves Random Media 10 (2000) 273–280. [9] S.H. Luo, L. Tsang, C.H. Chan, Application of the finite element method to Monte Carlo simulations of scattering of waves by random rough surfaces with the periodic condition, J. Electromagn. Waves Appl. 5 (8) (1991) 835–855. [10] Z. Cai, D. Chen, A direct searching method for nonlinear optimal problems and its convergence proof (in Chinese) (submitted for publication). [11] D.L. Jaggard, X. Sun, Scattering from fractally corrugated surfaces, J. Opt. Soc. Amer. A 7 (6) (1990) 1131–1139. [12] Z. Li, Y. Jin, Numerical simulation of bistatic scattering from fractal rough surface in the finite element method, Sci. China Ser. E 44 (1) (2001) 12–18. [13] Y.Q. Jin, Z. Li, Bistatic scattering and transmitting through fractal rough dielectric surface using FBM/SAA method, J. Electromagn. Waves Appl. 16 (4) (2002) 551–572. [14] Y.Q. Jin, Z. Li, Reconstruction of roughness profile of fractal surfaces from scattering measurement at grazing incidence, J. Appl. Phys. 89 (3) (2001) 1922–1926.