Inverse scattering solutions for side-band signals

Technological University Dublin ARROW@TU Dublin Conference papers School of Electrical and Electronic Engineering 2009-06-01 Inverse Scattering Solutions for Side-Band Signals Jonathan Blackledge Technological University Dublin, [email protected] Timo Hamalainen University of Jyvaskyla, Finland, [email protected] Jyrki Joutsensalo University of Jyvaskyla, Finland, [email protected] Follow this and additional works at: https://arrow.tudublin.ie/engscheleart Part of the Signal Processing Commons Recommended Citation Inverse Scattering Solutions for Side-Band Signals, Blackledge, Jonathan; ISSC 2009, UCD, June 10-11th 2009 , University College Dublin, 2009 This Conference Paper is brought to you for free and open access by the School of Electrical and Electronic Engineering at ARROW@TU Dublin. It has been accepted for inclusion in Conference papers by an authorized administrator of ARROW@TU Dublin. For more information, please contact [email protected], [email protected]. This work is licensed under a Creative Commons Attribution-Noncommercial-Share Alike 4.0 License ISSC 2009, UCD, June 10-11th Inverse Scattering Solutions for Side-Band Signals Jonathan Blackledge∗ , Timo Hämäläinen† and Jyrki Joutsensalo‡ SFI Stokes Professor of DSP Department of Mathematical Information Technology Dublin Institute of Technology University of Jyväskylä Ireland Finland E-mail: ∗ [email protected] † [email protected], ‡ [email protected] Abstract — When a signal is recorded that has been physically generated by some scattering process (the interaction of electromagnetic waves with an inhomogeneous dielectric, for example), the ‘standard model’ for the signal (i.e. information content convolved with a characteristic Impulse Response Function) is usually based on a single scattering approximation. An additive noise term is introduced into the model to take into account a range of non-deterministic factors including multiple scattering that, along with electronic noise and other background noise sources, is assumed to be relatively weak. Thus, the standard model is based on a ‘weak field condition’ and the inverse scattering problem is often reduced to the deconvolution of a signal in the presence of additive noise. Attempts at solving the exact inverse scattering problem which take into account multiple scattering effects often prove to be intractable, particularly with regard to the goal of implementing algorithms that are computationally stable and/or compatible with standard signal analysis methods and Digital Signal Processing ‘toolkits’. This paper provides an approach to solving the multiple scattering problem for narrow side-band systems (typically, electromagnetic signal processing systems) that is compounded in the introduction of a single extra term to the standard model. Keywords — Inverse Scattering, Multiple Scattering, Side-band Systems, Signal Reconstruction, Synthetic Aperture Radar I Introduction We consider a one-dimensional electromagnetic wave equation of the form (the inhomogeneous Helmholtz equation)  2  ∂ + k 2 u(x, k) = −k 2 γ(x)u(x, k) ∂x2 (1) where u is the electric wavefield, k is the wavenumber and γ = ǫr − 1, ǫr > 1 being the relative permittivity. The function ǫr is taken to be real whereas the wavefield u is taken to be a complex function (for generality) with variations as a function of x and k in both amplitude and phase. On the basis of this model, the forward scattering problem is defined as: Given γ obtain a solution for u. The inverse scattering problem is: Given u derive a solution for γ, e.g. [1], [2] and [3]. Let u(x, k), x ∈ (−∞, ∞) be given by the sum of an incident wavefield ui (x, k) and a scattered wavefield us (x, k), ui being given by a solution of  2  ∂ 2 + k ui (x, k) = 0. (2) ∂x2 Thus, with u = ui + us ,   2 ∂ 2 + k us (x, k) = −k 2 γ(x)[ui (x, k)+us (x, k)] ∂x2 where it is assumed that γ is of compact support, i.e. γ(x)∃∀x ∈ [−X, X]. The conventional approach to solving the inverse problem is as follows: (i) derive a solution for the scattered field us ; (ii) use this solution to obtain a solution for γ. The solution for us is [4], [5] us (x, k) = k 2 g(| x |, k) ⊗ γ(x)[ui (x, k) + us (x, k)] where ⊗ denotes the convolution integral and g is the Green’s function given by g(| x |, k) = i exp(ik | x |). 2k This equation requires an iterative procedure to be adopted in order determine us . However, under the weak field condition when kus k << kui k we consider the solution us (x, k) = k 2 g(| x |, k) ⊗ γ(x)ui (x, k). This is known as the Born approximation and has an asymptotic solution (for x → ∞ and with ui (x, k) = exp(−ikx)) given by [6] ik exp(ikx) us (x, k) = 2 Z∞ γ(x) exp(−2ikx)dx. −∞ Thus, in the far field, the relationship between the scattering function and the scattered field is characterised by the Fourier transform, a relationship that is specific to the weak field condition. The scattering amplitude spectrum As (k) (the ‘reflection coefficient’) can be defined as (with k → k/2) As (k) = where and γ e(k) = Z∞ ik γ e(k) 4 γ(x) exp(−ikx)dx −∞ 1 γ(x) = 2π Z∞ −∞ γ e(k) exp(ikx)dk The purpose of this paper is to explore a solution for γ and thereby us when kus k ∼ kui k that is based on the result, from equation (1):   u∗ 1 ∂2 γ(x) = − 1 + 2 2 us (x, k) | u(x, k) |2 k ∂x II Inverse Solutions Given equation (1), we can write   ∂2 1 k2 + u(x, k) γ(x) = − 2 k u(x, k) ∂x2 which, noting that γ = ǫr − 1 can be written in the form ǫr (x) = − k2 ∂2 u∗ (x, k) u(x, k) 2 | u(x, k) | ∂x2 and since u = ui + us we can write u∗i (x, k) + u∗s (x, k) k 2 | ui (x, k) + us (x, k) |2   ∂2 × k2 + us (x, k). ∂x2 γ(x) = − Note that when kus k << kui k,   2 ∂ 2 us (x, k) = −k 2 γ(x)ui (x, k), + k ∂x2 the equivalent expression for γ(x) being   u∗i (x, k) ∂2 2 γ(x) = − 2 k + 2 us (x, k) k | ui (x, k) |2 ∂x   2 ∗ ∂ u (x, k) k2 + us (x, k) =− i 2 k ∂x2 given that ui (x, k) = exp(±ikx) is a solution of equation (2). In comparison with the solution under the Born approximation, this solution for γ includes | u(x, k) |−2 which describes the (inverse square) amplitude envelop of the sum of the incident and scattered fields, i.e. with u(x, k) = Au (x, k) exp[iθu (x, k)], it is clear that | u(x, k) |2 = [Au (x, k)]2 . a) Solution for Side-band Signals If u is taken to be a narrow side-band signal that is dominated by a high frequency and carrier wave characterised by k0 , which, in turn, is determined by a narrow side-band incident wavefield ui , then we can consider the condition where | u(x, k0 ) |2 ∼ 1∀x (in general, a constant). This condition is equivalent to considering the case where u is a phase only field u(x, k0 ) = exp[iθu (x, k0 )] and is imposed in respect of the fact that, for narrow side-band signals, the contribution of | u |−2 to the reconstruction of γ from u is not significant when compared to (u∗i + u∗s )(k02 + ∂x2 )us . However, strictly speaking, the condition being imposed implies that, with ui = exp(±ik0 x), u∗i us + (u∗i us )∗ + | us |2 = 0. (3) The principal difference between the reconstruction of γ using the Born approximation and the inverse solution considered here is compounded in the addition of a single term, i.e. for a weak field where kus k << kui k   1 ∂2 ∗ γ(x) = −ui (x, k0 ) 1 + 2 2 us (x, k0 ) (4) k0 ∂x and for a strong field where kus k ∼ kui k   1 ∂2 γ(x) = −u∗i (x, k0 ) 1 + 2 2 us (x, k0 ) k0 ∂x   1 ∂2 −u∗s (x, k0 ) 1 + 2 2 us (x, k0 ). (5) k0 ∂x Note that under the strict definition of a phase only field which leads to equation (3), for kus k ∼ kui k γ(x) = ui (x, k0 )u∗s (x, k0 ) − [ui (x, k0 ) + us (x, k0 )]∗ ∂ 2 us (x, k0 ) k02 ∂x2 = −ui (x, k0 )u∗s (x, k0 ), k0 → ∞ In this paper, we relax this condition and utilze equations (4) and (5) given that | ui + us |2 ∼ 1. This is the ‘key’ to the results that follow. b) Numerical Simulation We consider a numerical simulation based on computing the scattered field using a forward differencing approach to equation (1). For a fixed wavenumber k0 (which defines a continuous wave mode), the vector un computed over a uniformly sampled discrete array composed of N samples, each of length ∆, is given by un+1 − 2un + un−1 +k02 un = −k02 γn un , n ∈ [1, N ] ∆2 which has the solution un = un+1 + un−1 un+1 + un−1 . = 2 2 2 − ∆ k0 (1 + γn ) 2 − ∆2 k02 ǫr,n (6) A solution to equation (6) can be obtained that is based on the following iteration: um+1 = n m um n+1 + un−1 2 − ∆2 k02 ǫr,n respectively where ⊗ now denotes the discrete convolution sum. Figures 1-2 show comparisons between the numerical results given by equations (8) and (9), illustrating the superiority of the inverse scattering solution considered over that given by the weak field solution for the case when k0 ∆ = 0.01. These results are based on applying the initial solution u1n = sin[2πip(n − 1)/(N − 1)] and Hilbert transforming the output arrays before computation of equations (8) and (9) through use of the MATLAB function hilbert. The profile of the reconstruction for ǫr,n based on equation (9) is preserved inclusive of characteristic ‘ringing’ due to the discontinuities associated with the model for ǫr,n . In comparison, the weak field solution given by equation (8) has a relatively narrow dynamic range and provides a poor reconstruction particularly with regard to resolving the discontinuities associated with ǫr,n . (7) where u1n is taken to be a discrete representation of the (unit amplitude) incident field which we define with a fixed number of periods p over n ∈ [1, N ], i.e.   2πip(n − 1) . u1n = exp (N − 1) A necessery condition that must be applied to equation (7) is that s 2 ∆k0 < kǫr,n k∞ where k • k∞ defines the ‘uniform norm’ and k0 = 2πp . N Fig. 1: Comparison between the weak and strong field inverse scattering solutions for the case when N = 10000, M = 100, ∆k0 = 0.01 with p = 50. From top to bottom: Relative permittivity function model ǫr,n ; real part of wavefield computed via equation (7); inverse solution (real part) computed using equation (8); inverse scattering solution (real part) computed using equation (9). Since a solution to the equation (1) is not conditional on the amplitude of the wavefield (i.e. the incident field can be A exp(±ik0 x) where A is an arbitrary value), the amplitude of the array um n, after M iterations is normalised on output, i.e. M M uM n → un /kun k∞ . For the discrete case, equations (4) and (5) transform to (for M iterations used to compute the scattered field) ǫr,n = 1 − u∗i,n [us,n + (∆k0 )−2 us,n ⊗ (1, −2, 1)] 1 −2 M = 1−(u1n )∗ [(uM (un −u1n )⊗(1, −2, 1)] n −un )+(∆k0 ) (8) and ǫr,n = 1 − u∗n [us,n + (∆k0 ) −2 us,n ⊗ (1, −2, 1)] ∗ M 1 −2 M = 1−(uM (un −u1n )⊗(1, −2, 1)] n ) [(un −un )+(∆k0 ) (9) Fig. 2: Comparison between the weak and strong field solutions for the case when N = 10000, M = 100, ∆k0 = 0.01 with p = 100. The descriptions of each plot follow those as given in Figure 1. III Solution for Narrow Side-Band Pulse-Echo Systems The inverse solutions derived and demonstrated numerically in the previous section have little practical value with regard to engineering a system designed to recover ǫr (x) from information on the scattered field measured ‘outside’ the scatterer, i.e. when | x |>| X |. This is because the solutions considered require that the scattered field is known ∀x including ∀x ∈ [−X, X]. Practically applicable systems typically measure the scattered field at a fixed point x in the far field (i.e. as x → ∞) by recording the spectrum of us over a range of values of k. In practice, this requires the application of pulse-echo mode methods. Pulse-echo methods involve the emission of a pulse and a recording of the back-scattered wavefield (echo). This approach is consistent with the physical nature a system if: (i) the scattering function is a one-dimensional function; (ii) the incident wavefield is a ‘pencil-line beam’. However, since all physical systems are intrinsically threedimensional, this model is idealised. Nevertheless, a variety of electromagnetic information and imaging systems can be ‘cast’ in terms of problems involving layered materials (e.g. the response of light, radio and microwaves to layered dielectric materials including the propagation of electromagnetic waves along transmission lines such as an optical fiber). For a narrow side-band system with carrier frequency k0 >>| k |, the inverse scattering solution is given by (for | Au |∼ 1)   ∂2 u∗ (x, k) 2 k0 + us (x, k) γ(x) = − k02 ∂x2 which, for a weak field where kus k << kui k reduces to   ∂2 u∗ (x, k) k02 + 2 us (x, k). γ(x) = − i 2 k0 ∂x These results apply for all values of x ∈ (−∞, ∞). As discussed in Section I, under the Born approximation, for x → ∞, a mapping is obtained between the scattered field us and γ that is based on the Fourier transform γ e(k) of γ(x), i.e. us (x, k) = ik0 exp(ikx)e γ (k), | k |<< k0 2 + u∗s (x, k0 )]  1 ∂2 1+ 2 2 k0 ∂x  us (x, k0 ) giving −[e u∗i (k, k0 ) + γ e(k) = u e∗s (k, k0 )]  | k |2 ⊙ 1− 2 k0  u es (k, k0 ) = −[e u∗i (k, k0 ) + u e∗s (k, k0 )] ⊙ u es (k), | k |<< k0 where γ e, u e∗i , u e∗s and u es are the Fourier transforms ∗ ∗ of γ, ui , us and us respectively and ⊙ denotes the correlation integral. Noting that, for ui (x, k0 ) = exp(−ik0 x), u e∗i (k, k0 ) = 2πδ(k0 − k) we have γ e(k) = 2πe us (k − k0 ) − u e∗s (k) ⊙ u es (k), | k |<< k0 or (ignoring scaing by 2π) u es (k) = pe(k)e γ (k + k0 ) + pe(k)[e u∗s (k + k0 ) ⊙ u es (k + k0 )], ∀k where pe(k) is some lowpass filter with a bandwidth significantly less than k0 . Thus, upon takling the inverse Fourier transform, we obtain an expression for the (demodulated) scattered field us given by us (x, k0 ) = p(x) ⊗ [γ(x) exp(−ik0 x)] +p(x) ⊗ [| us (x, k0 ) |2 exp(−ik0 x)] where p(x) is the inverse Fourier transform of pe(k). Here, p(x) described the bandlimited pulse that is incident on a layered dielectric described by γ(x). The first term is a description for the weak scattered field and the second term describes the effects of multiple scattering. In terms of the ‘standard model’ for a stationary signal s(t) as a function of time t where s(t) = p(t) ⊗ f (t) + n(t) (10) it is now clear that the Impulse Response Function (IRF) f is given by f (t) = γ(t) exp(−iω0 t) and the noise n(t) is given by n(t) = p(t) ⊗ [| s(t) |2 exp(−iω0 t)] which is a solution of u∗ (x, k) γ(x) = − i 2 k0 −[u∗i (x, k0 )  k02 ∂2 + 2 ∂x  us (x, k). This result suggests taking the Fourier transform of the equation γ(x) = where ω0 is the angular (carrier) frequency. Note that in practice, n(t) will include additional background noise and in this sense, we have extract a component of the noise term that is attributed to multiple scattering effects, albeit under the conditions that: (i) | ui + us |−2 ∼ 1; (ii) | k |<< k0 . IV Linear FM Signal Processing A linear FM ‘chirp’ signal, which is taken to be of compact support t ∈ [−T /2, T /2], is given by (in complex form and for a unit amplitude) p(t) = exp(iαt2 ), | t |≤ T 2 where α is a constant which defines the ‘chirp rate’ and T is the length of the ‘pulse’. The phase of this signal is given by αt2 (i.e. it has a quadratic phase factor) and its instantaneous frequency is therefore given by 2αt which varies linearly with time t. The frequency modulations are linear, the signal therefore being referred to as a ‘linear’ FM pulse. a) Inversion under the Weak Field Condition Consider a signal generated by the reflection of a linear FM pulse. Application of the weak field condition yields the following single scattering model: s(t) = exp(iαt2 ) ⊗ f (t), | t |≤ T . 2 Using a pulse of this type provides a simple but effective way of recovering the IRF by correlating the signal with the complex conjugate of the pulse to give the reconstruction T fˆ(t) = exp(−iαt2 ) ⊙ exp(iαt2 ) ⊗ f (t), | t |≤ 2 = T exp(−iαt2 )sinc(αT t) ⊗ f (t) ≃ T sinc(αT t) ⊗ f (t), T >> 1. This simplification, under the condition that T >> 1, is not only practically applicable, but systems desirable, especially with regard to electromagnetic pulse generation1 and allows the result for fˆ to be easily analysed in Fourier space. Using the convolution theorem we can write (ignoring scaling by π/α) ( F (ω), | ω |≤ αT ; F̂ (ω) = 0, | ω |> αT. which describes fˆ as being a band-limited version of f (given that f is not a band-limited function) where the bandwidth is determined by αT . b) Inversion Under the Strong Field Condition For a side-band FM system, the results discussed above imply that, under the single scattering approximation s(t) = p(t) ⊗ f (t) 1 A ‘long pulse’ provides greater signal energy and thus, the recorded data has a higher Signal-to-Noise Ratio. with inverse scattering estimate fˆ(t) = p∗ (t) ⊙ s(t) where, for a side-band system p(t) = exp(iαt2 ). Here, s(t) is taken to be the (complex) analytic signal obtained by taking the Hilbert transform of the real signal Re[s(t)]. For side-band systems this is typically undertaken by quadrature detection during demodulation from side-band to base-band [5]. The strong field inverse scattering solution for the signal, given by equation (10), implies that for the same (side-band) system fˆ(t) = T sinc(αT t) ⊗ f (t) +T sinc(αT t) ⊗ [| s(t) |2 exp(−iω0 t)]. V (11) Applications to Synthetic Aperture Radar Synthetic Aperture Radar (SAR) is a side-band pulse-echo system which utilizes the response of a scatterer as it passes through the radar beam to synthesize the lateral (azimuth) resolution [8], [9]. This allows relatively high resolution images to be obtained at a long range. The antenna emits a pulse of microwave radiation toward the ground and the return signal is recorded at fixed time intervals along the flight path. By studying the response of a SAR to a point scatterer in range x, and then in azimuth y, the Point Spread Function of the system is established which is given by [7] p(x, y) = sinc(αLx)sinc(βk0 y) where L is the pulse length and β is the angle of divergence of the radar beam. Thus, the (postprocessed) SAR image data fˆ(x, y) is given by the convolution of the object function f (x, y) (e.g. the ‘ground truth’) with the appropriate point spread function, i.e. fˆ(x, y) = p(x, y) ⊗2 f (x, y) where ⊗2 denotes the two-dimensional convolution integral and fˆ is taken to be complex data generated by f . A SAR image is usually generated by displaying the amplitude modulations of the data, i.e. ISAR (x, y) =| fˆ(x, y) | . The object function describes the imaged properties of the ground or sea surface. A conventional model for this function is the point scattering model that is typically used in SAR simulations [10]. This is where the object function is taken to be a distribution of point scatterers each of which reflects a replica of the emitted pulse and responds identically in azimuth. Here, nothing is said about the true physical nature of the been to work directly with the Helmholtz equation to produce an expression for the scattering function. The example numerical simulations given in Section II(b) provide evidence of the expected superiority of this solution over the weak scattering solution. However, it should be noted that this result is based on the application of the narrow side-band condition which is required in order to set | ui (x, k) + us (x, k) |2 ∼ 1. In this sense, the solution developed is not based on an entirely ‘exact’ inverse scattering solution, but rather a modified version of a solution tailored for applicability to narrow side-band systems. Such systems are 2 ˆ f (x, y) = p(x, y)⊗2 [f (x, y)+ | s(x, y) | exp(−ik0 x)] consistent with the applications of electromagnetic (12) signal processing. where the second term (on the right hand side) is A further extension of the inverse scattering aptaken to be the multiple scattering term2 . Figure proach considered here that is independent of the 3 shows the effect of applying a Gaussian lowpass condition | ui (x, k) + us (x, k) |−2 ∼ 1 can be develfilter to the complex data fˆ(x, y) before generation oped using the expansion of the image | fˆ(x, y) | using data obtained from the Sandia National Laboratories SAR database | ui (x, k)+us (x, k) |−2 = 1−ui u∗s −u∗i us − | us |2 +... [11]. Filtering the complex data is undertaken in order to attempt to suppress the term p(x, y) ⊗2 [| details of which lie beyond the scope of this paper. s(x, y) |2 exp(−ik0 x)]. The more usual practice is References to filter the amplitude image | fˆ |. However, within the context of equation (12), applying a lowpass [1] Z. S. Agranovich and V. A. Marchenko, ”The filter prior to computing an amplitude image of Inverse Problem of Scattering Theory”, Gorthe complex data reduces the influence of the cross don and Breach, 1963. terms inherent in | fˆ(x, y) | and it is this effect that [2] K. Chandan and P. C. Sabatier, ”Inverse reduces the speckle as illustrated in Figure 3. Problems in Quantum Scattering Theory”, Springer, 1977. object function in terms of material (dielectric) properties. Moreover, this standard convolution model for the data is based on application weak field approximation where multiple scattering effects (which generate ‘speckle’) are taken to be part of an additive noise function [7]. SAR images are coherent images (i.e. based on complex data containing magnitude and phase information) and consequently, contain noise that is characteristic of speckle patterns (coherent noise). Based on equation (11), we consider a model where [3] M. Cheney, J. H. Rose and B. DeFacio ”Three Dimensional Inverse Scattering”, Differential Equations and Mathematical Physics, Springer, 46-54, 1987. Fig. 3: Example of an original SAR image (left) and the same image after lowpass filtering the complex data (right). In both cases, the images have been histogram equalized utilizing the MATLAB function histeq. VI Conclusion The extension of the standard (stationary) model s(t) = p(t) ⊗ f (t) to include multiple scattering effects that are compounded in a single additional term provides a practical method of processing signals that is compatible with a standard signal processing ‘toolkit’. The inverse scattering problem is usually formulated by first solving the forward scattering problem. Once the relationship between the scattered field and the scattering function has been established by solving the Helmholtz equation, an inverse solution is then attempted. The approach taken in this paper has 2 After ing. application of conventional SAR signal process- [4] P. M. Morse and H. Feshbach, ”Methods of Theoretical Physics”, McGraw-Hill, 1953. [5] J. M. Blackledge, “Digital Signal Processing”, Horwood Publishing, 2006. [6] F, J. W. Olver, “Asymptotics and Special Functions”, Academic Press, 1974. [7] J. M. Blackledge, “Digital Image Processing”, Horwood Publishing, 2006. [8] R. O. Harger, “Synthetic Aperture Radar Systems”, Academic Press, 1970. [9] J. J. Kovaly, Artech, 1976. “Synthetic Aperture Radar”, [10] R, L. Mitchell, “Radar Signal Simulation”, MARK Resources, 1985. [11] http://www.sandia.gov/RADAR/sardata.html