Multidimensional wavefront estimation from differential ... - IEEE Xplore
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 39, NO. 3, MARCH 2001
655
Multidimensional Wavefront Estimation from Differential Delays Nicola Bienati and Umberto Spagnolini, Member, IEEE
Abstract—The propagating wavefield measured by complex sources-receivers arrangements can be characterized by a delayed waveform in the n-dimensional (n-D) space of measurements. In this paper, we propose a method for the nonparametric estimation of the time of delay function (also referred to as wavefront) that do not require the knowledge of the waveform. The n-D wavefront is estimated by using a least squares (LS) approach that integrates the local estimates of the differential times of delay of the waveform between neighboring sensors. The method shows a lower sensitivity to errors in local estimates when the overall multidimensionality (and redundancy) of the data is exploited. Examples considered here are derived from reflection seismology applications. Index Terms—Array signal processing, delay estimation, inverse problems, least squares methods, multidimensional signal processing, phase estimation.
I. INTRODUCTION
I
N REMOTE sensing, the features of the propagating medium can be investigated by measuring the backscattered wavefield (or simply the reflections) for various arrangements of sources and receivers. For acquisition efficiency (and data redundancy), the wavefield sampling is carried out by an array of receivers in various configurations. The measurements are 1)-D data volume where thus characterized by the ( spatial sampling (of sources and/or receivers) are located by . The wavefront originating from the n-D set each target is characterized by a delay function (or traveltime) in . Since the propagating medium is rather complex, the distortion with respect to the nominal wavefront is too difficult to be described by using a parametric method except for those distortions that leave the wavefront unchanged [1]. In this paper, we propose to estimate the delay function by using a nonparametric method (i.e., all the delays in are treated as free parameters). The only assumptions are that is continuous and single valued, the waveform is slowly changing across the sampling array, and the wavefield must be properly sampled compared to its minimum wavelength. For the wavefront estimation, many parametric approaches have been proposed in the literature (see e.g., [2] and references therein), while nonparametric strategies are known to be strongly nonlinear and unreliable. Routinely used nonparametric methods are based on simple ideas that hardly cope
Manuscript received June 1, 1999; revised March 21, 2000. The authors are with the Dipartimento di Elettronica e Informazione, Politecnico di Milano, I-20133 Milano, Italy (e-mail: [email protected]; [email protected]). Publisher Item Identifier S 0196-2892(01)02084-8.
with the aforementioned assumptions. The delays can be ) as reference estimated by taking one measurement ( and then by estimating all the differential delays of with respect to for . Since the waveform can change across the array, the differential delay estimates become are too far from less reliable when the measurement points . Alternatively, the differential delays can be estimated for neighboring points (or equivalently for those sensors that have similar waveform) or along a specified path and then all these . However, if estimates are collected into a unique solution the differential estimates are not constrained by a priori search, there is no guarantee that we will have a continuous and slowly changing wavefront. The main difficulties to constrain the wavefront continuity arise when dealing with multidimensional can be biased by the choice of the path arrays. The delay within the n-D array as the delay function along the sliced path is performed by optimizing the estimate locally along the path and not globally for the overall data. The key idea underlying the least squares (LS) wavefront estimation algorithm proposed here finds its root in the analogy between the two-dimensional (2-D) wavefront estimation and the 2-D phase unwrapping. Phase unwrapping can be considered as a wavefront estimation for monochromatic signals [3]. Similar in n-D can be estito 2-D phase unwrapping, the delay mated in two steps [4], [5]. First, the differential delays between neighboring sensors are estimated by using the cross-correlation method (in 2-D phase unwrapping this is referred to as instantaneous frequency estimation), then all of these differential deby implicitly constraining lays are integrated to estimate its smoothness and continuity. The main advantage is that each step can be optimized separately for the specific application. The methods for differential delay estimation are not covered in this paper as being extensively studied in the past, the reader is referred to [6], [7] and references therein. However, if differential delays are estimated from cross correlation of signals for couples of neighboring sensors, the errors may occur in picking the lag where the cross correlation peaks. These errors, also referred to as mispicks, are caused by low SNR and limited waveform bandwidth (e.g., when dealing with narrowband waveforms, the cross correlation becomes quasi-periodic and the delay estimate can skip one or more cycles). The way to is by considering compensate the effects of these errors on the wavefront estimation from differential delay estimation as a global (weighted or unweighted) LS problem. The research has been motivated by the need to estimate the wavefront or wavefront distortion (if estimation is performed after correction with nominal shape parameters) and therefore, some of the applications are described at first (Section II). In
0196–2892/01$10.00 © 2001 IEEE
656
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 39, NO. 3, MARCH 2001
Section III, the LS algorithm for an n-D array is proposed in a common framework for regularly and irregularly spaced sensors. The LS algorithm for 2-D array is discussed in Section IV. Since the wavefront estimation depends on the differential delays, the influence of mispicks on the overall solution has been evaluated. The analysis in Section V shows that errors due to mispicks spatially decrease faster when the dimensions of the array (or equivalently the redundancy of differential delay measurements) are increased. In addition, since the mispicks are either Gaussian and correlated or non-Gaussian, the weighted LS or the iteratively reweighted LS approach are more suitable in real applications. Performances have been evaluated in Section VI by taking a reflection seismology application. repreRemark on data dimensions: Measurements sents an ( 1)-D data volume, the n-D arrangement of sources and receivers is simply indicated as the n-D array of sensors (as for the reciprocity sources and receivers can be interchanged). When sources and receivers are placed everywhere in the 6. three-dimensional (3-D) space, the array of sensors is 4 (or the Sources and receivers arranged in planes lead to data volume can be, at most, five-dimensional [5-D]), while if 2 and the data sources and receivers are along lines, then volume is 3-D. The LS method is suitable to handle n-D array of sensors provided that the signals in the neighborhood of each sensor (either source or receiver) can be accessed. In seismic 4 [8]. acquisitions, arrays are planar and II. MOTIVATIONS AND APPLICATIONS In many remote sensing applications, the properties of the propagating medium, here described by a vector that contains the model parameters, are estimated from the minimization of delays the misfit between the estimated from for locations and the modeled according traveltimes (here is the to the model parameters matrix or vector transpose) (1) The minimization strategy and the stabilization of the objective , see e.g., function mostly depends on the nonlinearity of [9]. The method proposed in this paper has been motivated by the need to reliably estimate the traveltimes when the waveform is wideband (for velocity model estimation and static corrections computation in geophysics) or narrowband (for multidimensional phase unwrapping in digital elevation mapping by SAR interferometry). The nonparametric approach has been preferred to avoid any predefined relationship between the delay estimation and the optimization (1). In this section, we revise some applications that are based on the nonparametric estimates of the delays . In other cases not discussed here, there has been an attempt to avoid the explicit delays estimation in the minimization (1). A. Statics and Traveltime Tomography In reflection seismology, the near-surface velocity anomalies are known to distort the propagating wavefield. Such distortions (statics) contribute directly to the misfocusing of deep reflec-
tors in migrated gathers. In order to correct for these effects, it is useful to limit the analysis to kinematics, thus assuming that statics can be compensated by recovering only the differential delays between neighboring traces. These time shifts can be estimated simply by cross-correlating the traces after the compensation of the nominal (or undistorted) wavefront (this processing step is usually referred to as normal moveout correction). The key to the simplicity of this approach is that all the delays are estimated independently, but this is also its main weakness, as independent estimates may lack any continuity. In order to perform (residual) static correction that minimizes the propagation of wrong delay estimates, it is necessary to exploit some of the properties of in (1). For example, the differential time de(surface consistency) in order lays can be fitted to a model to exploit data redundancy for noise attenuation [10]. Alternatively, the time shifts can be estimated by maximizing the energy of the spatial average (or stack) of the signal, thus working directly on seismic data without any model constraint [11]. Even if the latter approach leads to computationally expensive optimizations of a nonconvex objective function, it is known to cope with severe statics [12]. All of these methods exploit the redundancy implied by the vertical straight rays propagation in the near surface. This assumption is reasonable only when 1) the velocity anomaly is very shallow compared to the target-reflector depth, and 2) the velocity of shallow layers is smaller than for deeper layers. Residual static corrections may fail if the aforementioned model assumptions are violated [13]. This occurs when velocity anomalies do not give rise to a constant delay for all events in each trace and they are so severe as to induce raypath bending, or when near surface lateral velocity variations are longer than the cable and fall in the null space of the inverse problem. In all these cases, results can be improved by estimating explicitly a model for velocity anomalies by means of nonlinear traveltime [15], [16]. Furthertomography [14] or linearization of more, velocity model building cannot be avoided when the final goal is prestack depth migration of the data. Similarly to statics reflection tomography is also affected by ambiguity between velocity and depth [17]. The advantage of tomographic algorithms is that this ambiguity can be easily controlled by constraining the velocity model through any available a priori information or through any other physically meaningful constraint. In this way, the resulting model can always be made reasonable [18]. In tomographic inverse problems, model fitting (1) is carried out by comparing the estimated traveltimes with those of rays traced through the model. The traveltime estimation is usually performed by manual interpretation of prestack data, and it is largely time consuming, especially for large n-D data volume from actual acquisitions. Therefore, it would be highly desirable to have a tool for automatic traveltime estimation that is faster and possibly more accurate than “human interpretation.” B. Phase Unwrapping In synthetic aperture radar (SAR) interferometry, two (or more) backscattered fields that pertain to the same scene recorded (not necessarily at the same time) by two (or more) SAR antennas allow us to determine the position of the scattering points. The delays between different travelling paths can
BIENATI AND SPAGNOLINI: MULTIDIMENSIONAL WAVEFRONT ESTIMATION
be used to obtain the elevation map of the scene [19]. When the waveform is narrowband or monochromatic, the delays can be estimated only modulo the period of the center frequency (time dependency is (or equivalently modulo 2 ) from unnecessary and it can be dropped). In this case, the elevation map of a 2-D scene is obtained from the traveltimes only after phase unwrapping and re-scaling. 2-D phase unwrapping algorithms are based on the estimates of the differential phase (or their modifications); the phase difference for neighboring and this estimate is constrained pixels is measured modulo . The 2-D phase unwrapping to be within the interval algorithms are based on the integration of several local measurements of differential phases, each performed for the nearest neighbor of each measuring point [5]. Multiple interferograms ) phase give rise to the need of a multidimensional ( unwrapping strategy when there is the need to track the phase variations across the interferograms [20]. This problem has to deal with the simultaneous unwrapping of multiple 2-D planes, and it can be approached in a general framework as described below. Moreover, the LS approach can easily cope with irregularly sampled data. Recently, phase unwrapping on a sparse grid of stable targets has gained increasing interest (see, e.g., [21], [22]), and the use of the multidimensional approach improves considerably the reliability of the results [23]. More measurements of the same image can be carried out with the purpose of improving the SNR (i.e., the phase image shows no variation across the interferograms except for different noise realizations). In this case, the averaging of those solutions obtained by unwrapping each 2-D plane separately is equivalent to the solution that can be achieved by globally unwrapping the 3-D data volume with the constraint of having the same solution on each plane. This special case is not considered here. III. LS WAVEFRONT ESTIMATION After definition of the model (Section III-A), the LS algorithm is discussed for n-D array regardless of sensors’ positioning (Section III-B). The LS method for regularly spaced 2-D array is rederived in Section III-C by using a compact notation for structured matrixes. This is the same problem commonly solved in LS approach to 2-D phase unwrapping. Continuous-space analogy for 2-D [24] is instrumental to derive the properties of LS method for n-D array (Section IV). A. Data model 1)-D data volume Let the ( layed waveform
be described as a de-
(2) has been (approxiwhere the target event is any mately) windowed from the whole data volume, spatially and temporally structured noise (i.e., all the interfering influences only the events), and the amplitude variations for the target reflection (or for any SNR. The waveform compound reflection) need not be known and can also be slowly , with changing in space: . In other words, for any neighboring
657
point , the normalized correlation (or the coherence between the two waveforms)
(3)
1. The delay function is should be assumed to be a continuous and regular surface (up to the second order derivatives), and the continuity also implies that has to be single valued (without any triplications or faults). In the presence of diffractions, it may be necessary to backpropagate the wavefield (migration), estimate the delay , or in the migrated domain, and then demior , see e.g., [25]. grate the estimated delay B. LS algorithm Let arg
(4)
and estibe the differential delay between the signals in mated by using the cross-correlation method, the unknown travcan be estimated by solving, in the LS sense, eltime surface for the equality (5) measurements. The estimate within a predefinite subset of is thus given by the solution of the linear system (6) equations for unknowns (usually ), is the of delays, and is a vector consisting of the difvector of ferential delays between all the couples of traces useful for the differential delay estimation. is a sparse matrix, each row has only two nonzero elements that corresponds to the differential delay between the two sensors. The LS solution of the system (6) is known to be given by (7) and it can be computed efficiently by means of iterative methods that exploit the sparsity of (e.g., successive over-relaxation, preconditioned conjugate gradient, or QR decomposition; see [26] or [27] for a detailed discussion). The system (6) and its LS solution (7) represent the general formulation of the problem for n-D arrays, as the regularity of the array geometry is not required. Only the selection of the couples of sensors for the set of differential delay estimates is needed. The preferred choice is for those couple of traces where the waveform is slowly changing and the signal is mostly spatially-coherent, i.e., 1. In practice, the subset is lower than the sensors: ( 1)/2. overall number of couples for is in order. The A discussion on the structure of matrix can be carried out by sensors’ pairing to obtain the matrix considering a neighboring selection criteria. The subset of differential delays can be estimated by looking for the sensors pairs
658
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 39, NO. 3, MARCH 2001
and so that the geometrical distance , for . Depending on the array geometry, the number some value of measurements is linearly increasing with the number of sen, where is the number of neighbors to each sors: sensor. For regularly spaced sensors, the selection criteria for (the system surrounding sensors is trivial and it yields 2). In addition, when the ratio (6) is overdetermined only if between the number of constraints and the unknowns is constant within a large support (or in the whole support as for regularly spaced array), the solution is expected to have predictable properties. For an irregular acquisition geometry, the problem of determining the nearest neighbor of each source and/or receiver can be solved by using any of the known methods. However, here we preferred the use of the Delaunay triangulation algorithm [28] even if the sensors’ pairing algorithm is not critical as far as the selection criteria is reasonable. The structure of matrix depends only on the sensors ordering. The properties and the limitations of the LS approach can be evaluated only for some simple (but useful) cases with regular array geometry: the 2-D array and the n-D array. Most of the analytic properties of the LS approach can be more easily derived when considering the continuous-space wavefront (not sampled). This will be discussed in Section V.
is a small constant value that stabilizes the iterations. The . LS solution (7) represents the first iteration, i.e., It can be shown that this choice of the weighting matrix leads the IRLS to have performances similar to the L1-norm [29] and thus, IRLS turns out to be quite effective in reducing the influence of large mispicks in differential delays . IV. WAVEFRONT ESTIMATION FOR 2-D ARRAY For a uniformly spaced 2-D array, the differential delay for neighboring traces can be evaluated along the two axes here redenote the maferred to as and . Let R and trix of (unknown) delays, R the matrices of differential dethe lays for neighboring traces (or gradients), and shift matrix, i.e., where is the Kronecker 1 for 0, and 0 for 0). The matrix operdelta ( ators to compute the gradients from the matrix can be written and , where as is the unitary matrix. Let and be the differential delay matrix augmented by or zeros, respectively. The delays can be oba vector of tained from the following linear relationships on the gradients: (9)
C. Weighted and iteratively re-weighted LS The equality (5) can be weighted according to any predefinite criteria. For a positive definite weighting matrix (not necessarily diagonal), the minimization of the weighted leads to the following solution: norm . In addition to the weighting matrix , it is possible, when the Gaussian model holds true, to introduce the a priori information through the covariance matrix of traveltimes. Since these matrices are often unknown, feasible approximations are sometimes needed. The covariance matrix of can be chosen so that the dispersion of each measurement is proportional to the geometrical distance , or equivis diagonal with entries set as . This was alently adopted in the seismic applications discussed in Section II-A. with 0 1 can be used to Diagonal matrix pair the sensors according to any selection criteria based on the stationarity of either the wavefield or the SNR. In principle, at the beginning, all possible sensor pairs can be taken into ac( 1)/2. Then, after some iterations, terms are count and 0. In this latter case, downweighted or can be set to the corresponding equations are simply removed from system (6). can also be computed from The weighting matrix , and the LS solution (7) can be the residuals solved iteratively. This technique is also known as the iteratively be the estireweighted least squares (IRLS). In practice, let mated wavefront at the th iteration according to solution of the , weighted LS system for the next iteration depends on the the weight matrix , and it is diagonal overall residual
(8)
vec According to the equality vec for the vec( ) operator that stacks all columns of a matrix into one single column vector and some of the properties of Kronecker products [30], the system of (9) can be vec , rewritten similar to system (6) by defining vec vec and . , and the equaThe number of unknowns is 2 . For tions are matrix regularly sampled 2-D array, the is structured and sparse similar to Hunt [31]. The normal equations can be rewritten as (10) From the property of the Kronecker products, the normal equations [(10)] are equivalent to the discrete version of the Poisson’s equation comprehensive of the Neumann boundary (or ) conditions [3]. For instance, the matrix is the discrete version of the second order derivative along (or ). is identically zero, as it is associOne eigenvalue of ated to the eigenvector corresponding to constant shifts (this is expected as the LS wavefront estimation is based on the differential delays). For 2-D array, the error analysis becomes easy if the estimation problem is reconsidered in continuous space. Assume that we have estimated the gradient field (where and are the unit norm vectors for the and axes, respectively) of the delay function within the region of support . The wavefront can be estimated from the and , the minimization of the misfit between
BIENATI AND SPAGNOLINI: MULTIDIMENSIONAL WAVEFRONT ESTIMATION
gradient of (here ). The solution of the variational problem is provided by the integration of the Poisson’s equation (11) and are the Laplace and the divergence opwhere erators, respectively. The numerical solution of (11) with Neumann boundary conditions obtained by approximating partial derivatives with finite differences leads to (10). When (11) is used for the wavefront estimation, it has a nonempty data null space. This null space contains those components of the differential delays that cannot be predicted by the model, and it represents a residual misfit after the minimum has been achieved. This is due to the choice of describing the gradients field by using the scalar function only. According to the Helmoltz’s theorem, any vector can be represented by the superposition of a scalar field and a vector potential potential (12) cannot give rise to where the continuous surface itself , which is accounted for by the rotational component of the vector potential. This component must necessarily be part of 0 for any . the data null space as The rotational component is due to those values in the estimated gradients that conflict with the continuity of the solution caused either by mispicks and/or by not complying with the continuity hypothesis (e.g., when the surface is faulted or represents a linear transformamultivalued). Since matrix tion, the rotational component is also given by the projection of the measurements onto the null space of (see [3] for further details). When dealing with irregularly spaced array, the evaluation of the projection matrix is the only way to establish the sources of error. In any case, similar to 2-D phase unwrapping, the rotational component field represents a useful indicator of the quality of the cross-correlation maxima picking and of the match of the delay surface with the continuity hypothesis. V. ADVANTAGES OF WAVEFRONT ESTIMATION MULTIDIMENSIONAL ARRAYS
IN
In principle, n-D array can be sliced to extract several 2-D arrays and have simpler wavefront estimation algorithms. Here we demonstrate, by simple arguments, that this approach is far from being an optimum one except for trivial cases. A multidimensional approach for n-D array is expected to have a better constraint as, in the simpler case of regularly sampled arrays, the equations (one delays are obtained by solving a system of unknowns. Therefore, equation for each dimension) for the 3), a higher redundancy of mulfor ( 1)-D data volumes ( tidimensional algorithms guarantees better performances when compared to iterated estimation on lower dimensional subvolumes obtained by slicing the global volume. The advantages for n-D array can be evaluated in terms of error propagation by considering the continuous space case. Let contain only the estimated gradients one single incorrect measure due to a mispick (for instance in along one component only, say :
659
), the wavefront error caused by follows from the integration of the partial differential equation (13) , it is derived from (11) as for the true wavefront . On each dimension, the finite differences approximation splits the mispick into two neighboring terms of of opposite sign: a dipole. Since the driving term in (13) is a dipole, the solution is the combination of the responses of two pulses having opposite signs that decays faster than each pulse response. can be evaluated first by evaluating the The decaying of to one pulse: . response of for According to Gauss’s law, the amplitude of decays as the inverse of the surface of the [32]. For a dipole having unit width n-D sphere: spacing the decays is the combination of both terms (when is ) larger than the dipole width, in this case where even (14) odd or equivalently [by assuming that for radial]
the field (14) is
(15) In conclusion, the advantage of the multidimensional traveltime 2 is the lower spreading of those errors estimation for that lead to conflicting values in the gradient (also referred to as rotational sources in 2-D phase unwrapping [4]), these errors are mainly due to mispicks. This qualitative proof that mixes continuous and discrete doin (15) is given mains is far from being a rigorous one as consists of by the solution of (13) when the forcing term one pulse and it is measured far from the pulse-location. For the discrete-space the error behavior depends on the numerical solution of (7) that cannot be recursively computed as its support is not wedge shaped, the results are shown in Fig. 1 for 2-D and 3-D arrays. In both cases the mispick occurs on , the 3-D response is shown only on one plane as it is invariant by rotation. According to the relationship (15) the error (asymptotically) defor 3-D arrays. Moreover, also the cays as 1 for 2-D and 1 peak amplitude decreases as the number of spatial dimensions increases (from 1/4 for 2-D to 1/6 for 3-D arrays). When the errors are far from being Gaussian distributed (as for errors that arise from mispicks), the approximate L-1 norm approach or IRLS is advisable. VI. EXAMPLES AND PERFORMANCES FROM REFLECTION SEISMOLOGY In reflection seismology, the space-axes represent any suitable combination of space coordinates from the acquisition ge-
660
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 39, NO. 3, MARCH 2001
Fig. 1. Mispick response (x) for wavefront estimation in two-dimensional (2-D) array (left: the response has been magnified by a factor of 4) and in three-dimensional (3-D) array (right: the response has been magnified by a factor of 6).
Fig. 2.
Source-receiver geometry for 3-D land acquisiton (SEG/EAEG 3-D Overthrust B-3 acquisition subset).
ometry (source and receiver positions, midpoint and offset, in 2-D experiments; receiver positions in a common shot gather, offset-azimuth or any other suitable combination for 3-D acquisitions). Even if the LS wavefront estimation algorithm is general, here we tailor the method for 3-D and four-dimensional (4-D) data volumes. A short description of the 3-D acquisition geometry will be provided at first (for a detailed descriptions of source/receiver arrangements the reader is referred to [8]). The way how the source/receiver arrangements sample the wavefield depends on the specific acquisition templates [33]. denotes the coordinate associated to the In this section ) those associated to the receivers. For field sources and ( acquisitions the five-dimensional (5-D) data volume consists of sources and receivers, the overall 4-D array samples the wavefield into locations. The 3-D land acquisition geometry considered here has multiple cross-spreads (each cross-spread is a 3-D data volume as sources and receivers moves along orthogonal lines), as shown , in Fig. 2. The 4-D volume of the gathered data becomes and . The data considered for here
the examples proposed here are from the SEG/EAEG 3-D Overthrust B-3 acquisition subset (dataset available at [34], see [35] for a description of the experiment). According to Fig. 2, we selected eight consecutive source-receiver lines (cross-spreads) that correspond either to eight traveltime surfaces each obtained for a given or to an overall traveltime volume . The main feature of this acquisition is the nonuni325 m 50 m form spacing of receiver lines where for , 8. The sensors are uniformly spaced along the receiver lines (spacing is 50 m). The reflections of the target event arise from a nonflat reflector embedded in complex interfering backscattered wavefields (data and the estimated volume are in Fig. 3. This subset is realistic enough to be used here to compare the performances of the LS method for varying SNR, increasing dimensions, and sparse/irregular acquisitions. IRLS is also evaluated. The advantages of the wavefront estimation for multidimensional arrays is shown in Fig. 4, which quantifies, for the data 2 receiver lines in 350 m volume obtained from 400 m, the better immunity to noise when dealing with and
BIENATI AND SPAGNOLINI: MULTIDIMENSIONAL WAVEFRONT ESTIMATION
661
Fig. 3. Slices from the subset of SEG/EAGE three-dimensional (3-D) Overthrust B-3 data. 3-D data volume d(x ; x = 325m; y ; t) sliced along three axes [planes in gray-scale: d(x ; x = 325 m, y ; t = 0.85 s), d(x = 2000 m, x = 325 m, y ; t), and d(x ; x = 325 m, y = 900 m, t); white lines: signals in four receiver-locations] and the estimated delays (x ; x ; y ) for the first cross-spread (x = 325 m). Wavefront correction has been applied for visual display to flatten the nearly hyperbolic wavefront.
Fig. 4. Root mean square error (RMSE) of estimated delays versus noise percentage for 2-D array (or 3-D data volume) and 3-D array (or 4-D data volume). Experimental results after 50 runs (N = 2, N = 120, N = 31, N = 1).
4-D data volumes compared to 3-D. The comparison is given in terms of root mean squared error (RMSE) of estimated traveltime for increasing levels of (white Gaussian) noise superimposed to the data set. The four surfaces (estimated separately on each 3-D data volume) and the volume (estimated globally from the 4-D data volume) obtained in the absence of Gaussian noise have been used as references to evaluate the RMSE after the noise has been added. The standard deviation of the noise level was normalized to the root mean squared value of the signal as this can be more easily related to the way the noise level is
Fig. 5. Wavefront estimation (x ; x ; y ) for 3-D array (cross-spread acquisition). The 3-D delay surface is estimated by using the full 4-D data volume are plotted for varying x = {325 m, 625 m, 925 m, 1225 m} (view rotated by 90 with respect to Fig. 3).
inferred from experimental measurements. In Fig. 6 the same comparison between the strategies that deal with 3-D and 4-D data volumes is shown for 95% and 160% of noise. In any case, for noise below 120% the RMSE is approx. below the sampling interval (8 ms). The advantages of higher dimensions, mainly for higher noise level, are justified by the larger redundancy. In addition, the advantages of the IRLS approach for large noise is a consequence of the interaction between the noise and the cycle-skips for low SNR.
662
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 39, NO. 3, MARCH 2001
=
Fig. 6. Slices of data volume (same as Fig. 3) for the first cross-spread (x 325 m) and estimated wavefront for increasing additive noise level. Left column: the wavefront is estimated in 2-D array. Right column: the wavefront is estimated in 3-D array (the same slice as for 2-D array is shown in right column).
When the estimation is carried out through subvolumes slicing, the estimated wavefronts can be mistied as each estimate cannot be easily tied together with the others (in some cases, the slicing strategy can be arranged to have several common points among different slices). For example, in acquisition comprising multiple cross-spreads (see Fig. 2) each traveltime surface can be estimated from each cross-spread can be (2-D array) or the overall delay volume estimated from the whole data set. Once estimated, the volume can be sliced into surfaces if necessary, but each of these surfaces are naturally tied by multidimensional estimation and they are not floating as for traveltime estimation after slicing (where the relative shift of each solution is arbitrary). This is another reason to exploit the full dimensionality of the data set. Fig. 5 shows the four traveltime surfaces (for 325 m, 625 m, 925 m, 1225 m) obtained by the 4-D algorithm from the mentioned subset of SEG/EAEG 3-D Overthrust. Due to the slope of the target reflector the surfaces have a different shape but 4-D estimation has correctly tied the delay surfaces shows an increasing where the reflector is flat and delay when the receiver lines move up-dip. The same data set considered in Fig. 6 has been randomly decimated 2:3 and 2:5 in order to simulate an acquisition with sparse source-receivers randomly located. The delay function
has been estimated on a single decimated cross-spread (2-D array). Fig. 7 shows that even in absence of additive noise the reduction of the number of traces increases the error in the estimated wavefront due to a lower coherence (3) between the waveform across the array caused by a larger distance between neighboring sensors. Obviously the results worsen when white Gaussian noise is added to the data (Fig. 7, right). In order to keep the same redundancy for different decimations, the differential delays have been estimated on a number of pairs so that the ratio is approximately con(or equations) stant. The advantage of the IRLS method is shown in Fig. 8 for 4-D data volumes and increasing noise level. Here the first 3 receiver lines have been considered having sensors’ {325 m, 375 m, 625 m}, respectively. The spacing of greater distance between the second and the third line (250 m) implies that the slowly changing waveform hypothesis is not fully satisfied, and it increases the probability of cycle-skip occurrence in cross-correlation maxima picking. Even if the large array spacing makes the advantages of the multidimensional approach rather negligible, the wavefronts are tied consistently (e.g., Fig. 5). In addition, the interaction between high noise level and large errors (mainly due to cycle-skips in differential delay estimation) makes the IRLS method an attractive tool to handle low SNR data.
BIENATI AND SPAGNOLINI: MULTIDIMENSIONAL WAVEFRONT ESTIMATION
663
=
Fig. 7. Wavefront estimation in 2-D array for decimated data extracted from the first cross-spread (x 325 m). Left: wavefront estimation increasing decimation rate (2:3 and 2:5, noise free data). Right: wavefront estimation for noisy data (95% of Gaussian noise) and increasing decimation rate.
Fig. 8. RMSE of estimated delays for 3-D array (irregular spacing for x ) versus noise percentage for LS approach and IRLS. Experimental results after 50 runs ( 3, 120, 31, 1).
N = N =
N =
N =
VII. CONCLUSIONS This paper has shown that a least squares approach (or any modifications such as weighted LS or iteratively reweighted
LS) can be used for the nonparametric estimate of a delay (hyper)volume from the differential delays for neighboring sensors. Cross-correlation maxima picking is used as differential delays estimators. The errors due to peak selection ambiguity of cross-correlator estimator (mispick) have been evaluated and their propagation has been proved to be limited by the LS method that exploits the spatial continuity of the wavefront over the multidimensional space. The robustness against noise and mispicks improves with the dimensionality of the data set. The L1-norm approach is advisable when the mispicks are more likely (e.g., when the sensor’s spacing is too large and the waveform changes across the array). In the proposed method, the constraint that the differential estimates are only for neighboring sensors has been relaxed. Therefore, it can be extended to the differential delay estimates for any sensors’ configuration (and thus for irregularly spaced sensors) and any waveform (wideband or pulsed waveforms and narrowband or phase unwrapping problems). ACKNOWLEDGMENT The authors would like to thank F.Rocca for some stimulating discussions on the application of the multidimensional approach to geophysics.
664
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 39, NO. 3, MARCH 2001
REFERENCES [1] U. Spagnolini and S. Grion, “Shape parameter estimation of wavefronts with known waveform,” Signal Process., vol. 68, pp. 233–257, Aug. 1998. [2] D. H. Johnson and D. E. Dudgeon, Array Signal Processing: Concepts and Techniques. Englewood Cliffs, NJ: Prentice-Hall, 1993. [3] U. Spagnolini, “Non-parametric narrow-band wavefront estimation from wavefront gradients,” IEEE Trans. Signal Processing, vol. 47, pp. 3116–3121, Nov. 1999. , “2-D phase unwrapping and instantaneous frequency estimation,” [4] IEEE Trans. Geosci. Remote Sensing, vol. GRS-33, pp. 579–589, May 1995. [5] D. C. Ghiglia and M. D. Pritt, Two-Dimensional Phase Unwrapping. New York: Wiley, 1998. [6] IEEE Trans. Acoust., Speech, Signal Processing, vol. ASSP-29, Pt. II, June 1981. [7] G. C. Carter, “Coherence and time delay estimation,” Proc. IEEE, vol. 75, pp. 236–255, Feb. 1987. [8] D. G. Stone, Designing Seismic Surveys in Two and Three Dimensions. Tulsa, OK: Soc. Expl. Geophys., 1994. [9] J. A. Scales, P. Docherty, and A. Gersztenkorn, “Regularization of nonlinear inverse problem: Imaging the near-surface weathering layer,” Inv. Probl., vol. 6, pp. 115–131, 1990. [10] R. A. Wiggins, K. L. Larner, and +R. D. Wisecup, “Residual statics analysis as a general linear inverse problem,” Geophysics, vol. 41, no. 5, pp. 922–938, Oct. 1976. [11] J. Ronen and J. F. Claerbout, “Surface consistent residual statics estimation by stack-power maximization,” Geophysics, vol. 50, no. 12, pp. 2759–2767, Dec. 1985. [12] D. H. Rothman, “Nonlinear inversion, statistical mechanics, and residual statics estimation,” Geophysics, vol. 50, no. 12, pp. 2784–2796, Dec. 1985. [13] D. Palmer, Refraction Seismics. London, U.K.: Geophysical Press, 1986. [14] T. N. Bishop, K. P. Bube, R. T. Cutler, R. T. Langan, P. L. Love, J. R. Resnick, R. T. Shuey, D. A. Spindler, and H. W. Wyld, “Tomographic determination of velocity and depth in laterally varying media,” Geophysics, vol. 50, pp. 903–923, June 1985. [15] E. Loinger, “A linear model for velocity anomalies,” Geophys. Prospect., vol. 31, pp. 98–118, Feb. 1983. [16] F. Rocca and U. Spagnolini, “A simple interpretation of reflection tomography,” in Proc. Eur. Assoc. Exploration Geophysicists (EAEG) Meeting, vol. 1, Glasgow, U.K., 1995, pp. D2–D3. [17] C. Stork, “Singular value decomposition of the velocity-reflector depth tradeoff, Part 2: High-resolution analysis of a generic model,” Geophysics, vol. 57, pp. 933–943, July 1992. [18] C. Stork and R. W. Clayton, “Using constraints to address the instabilities of prestack velocity analysis,” Geophysics, vol. 57, pp. 404–419, Mar. 1992. [19] H. A. Zebker and R. M. Goldstein, “Topographic mapping from interferometric synthetic aperture radar observations,” J. Geophys. Res., vol. 91, pp. 4993–4999, 1986. [20] A. Ferretti, C. Prati, F. Rocca, and A. Monti Guarnieri, “Multi-baseline SAR interferometry for automatic DEM reconstruction,” in Proc. 3rd ERS Symp., Florence, Italy, 1997. [21] M. Costantini and P. Rosen, “A generalized phase unwrapping approach for sparse data,” in Proc. Int. Geosci. Remote Sensing Symp., Hamburg, Germany, 1999, pp. 267–269. [22] M. Eineder and J. Holzner, “Phase unwrapping of low coherence differential interferograms,” in Proc. Int. Geosci. Remote Sensing Symp., Hamburg, Germany, 1999, pp. 1727–1730.
[23] A. Ferretti, C. Prati, and F. Rocca, “Non-linear subsidence rate estimation using permanent scatterers in differential SAR interferometry,” IEEE Trans. Geosci. Remote Sensing, vol. 38, pp. 2202–2212, Sept. 2000. [24] H. Takajo and T. Takahashi, “Noniterative method for obtaining the exact solution for the normal equation in least-squares phase estimation from the phase difference,” J. Opt. Soc. Amer. A, vol. 5, pp. 1818–1827, 1988. [25] J. A. C. Jacobs, F. Delprat-Jannaud, A. Ehinger, and P. Lailly, “Sequential migration-aided reflection tomography: A tool for imaging complex structures,” in 62nd Annu. Int. Meeting SEG, Expanded Abstract, New Orleans, LA, 1992. [26] G. H. Golub and C. F. Van Loan, Matrix Computations, 2nd ed. London, U.K.: Johns Hopkins Univ. Press, 1989. [27] G. Meurant, Computer Solution of Large Linear Systems. Amsterdam, The Netherlands: Elsevier, 1999. [28] F. Preparata and M. Shamos, Computational Geometry—An Introduction. Berlin, Germany: Springer-Verlag, 1985. [29] J. A. Scales and A. Gersztenkorn, “Robust method in inverse theory,” Inv. Probl., vol. 4, pp. 1071–1091, 1988. [30] A. Graham, Kronecker Products and Matrix Calculus. New York: Wiley, 1981. [31] B. R. Hunt, “Matrix formulation of the reconstruction of phase values from phase differences,” J. Opt. Soc. Amer. A, vol. 69, no. 3, pp. 393–399, Mar. 1979. [32] R. Courant and T. Hilbert, Methods of Mathematical Physics. New York: Wiley. [33] G. J. O. Vermeer, Seismic Wavefield Sampling. Tulsa, OK: Soc. Expl. Geophys., 1990. [34] SEG/EAEG 3-D Modeling data access site: http://archiv.llnl.gov/SSD, . [35] F. Aminzadeh, N. Burkhard, P. Lailly, F. Rocca, and K. Wyatt, “Progress Rep. from the SEG/EAEG 3-D modeling commitee,”, vol. 13, 1994.
Nicola Bienati received the Dott. Ing. degree (cum laude) in telecommunications and the Ph.D. degree in electronics and telecommunications from the Politecnico di Milano, Milano, Italy, in 1995 and 2000, respectively. His main research interests are seismic reflection tomography, horizon picking in multidimensional data volumes, blind deconvolution, and blind sources separation.
Umberto Spagnolini (M’99) received the Dott. Ing. Elett. degree in telecommunications (cum laude) from the Politenico di Milano, Milano, Italy, in 1988. Since 1988, he has been with the Dipartimento di Elettronica e Informazione, Politenico di Milano, where he has held the position of Associate Professor of Digital Signal Processing since 1998. His primary research (and applications) focuses on array processing and wavefield interpolation (mobile communication and geophysics), inverse problems (ground penetrating radar) and parameter estimation (2-D phase unwrapping for SAR), and non-Gaussian EMI reduction. Dr. Spagnolini is a member of the SEG and EAGE and serves as an Associate Editor for IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING. He was awarded the Associazione Elettrotecnica ed Elettronica Italiana (AEI) Award and the Van Weelden Award of EAGE, both in 1991, and received the Best Paper Award from EAGE in 1998.
655
Multidimensional Wavefront Estimation from Differential Delays Nicola Bienati and Umberto Spagnolini, Member, IEEE
Abstract—The propagating wavefield measured by complex sources-receivers arrangements can be characterized by a delayed waveform in the n-dimensional (n-D) space of measurements. In this paper, we propose a method for the nonparametric estimation of the time of delay function (also referred to as wavefront) that do not require the knowledge of the waveform. The n-D wavefront is estimated by using a least squares (LS) approach that integrates the local estimates of the differential times of delay of the waveform between neighboring sensors. The method shows a lower sensitivity to errors in local estimates when the overall multidimensionality (and redundancy) of the data is exploited. Examples considered here are derived from reflection seismology applications. Index Terms—Array signal processing, delay estimation, inverse problems, least squares methods, multidimensional signal processing, phase estimation.
I. INTRODUCTION
I
N REMOTE sensing, the features of the propagating medium can be investigated by measuring the backscattered wavefield (or simply the reflections) for various arrangements of sources and receivers. For acquisition efficiency (and data redundancy), the wavefield sampling is carried out by an array of receivers in various configurations. The measurements are 1)-D data volume where thus characterized by the ( spatial sampling (of sources and/or receivers) are located by . The wavefront originating from the n-D set each target is characterized by a delay function (or traveltime) in . Since the propagating medium is rather complex, the distortion with respect to the nominal wavefront is too difficult to be described by using a parametric method except for those distortions that leave the wavefront unchanged [1]. In this paper, we propose to estimate the delay function by using a nonparametric method (i.e., all the delays in are treated as free parameters). The only assumptions are that is continuous and single valued, the waveform is slowly changing across the sampling array, and the wavefield must be properly sampled compared to its minimum wavelength. For the wavefront estimation, many parametric approaches have been proposed in the literature (see e.g., [2] and references therein), while nonparametric strategies are known to be strongly nonlinear and unreliable. Routinely used nonparametric methods are based on simple ideas that hardly cope
Manuscript received June 1, 1999; revised March 21, 2000. The authors are with the Dipartimento di Elettronica e Informazione, Politecnico di Milano, I-20133 Milano, Italy (e-mail: [email protected]; [email protected]). Publisher Item Identifier S 0196-2892(01)02084-8.
with the aforementioned assumptions. The delays can be ) as reference estimated by taking one measurement ( and then by estimating all the differential delays of with respect to for . Since the waveform can change across the array, the differential delay estimates become are too far from less reliable when the measurement points . Alternatively, the differential delays can be estimated for neighboring points (or equivalently for those sensors that have similar waveform) or along a specified path and then all these . However, if estimates are collected into a unique solution the differential estimates are not constrained by a priori search, there is no guarantee that we will have a continuous and slowly changing wavefront. The main difficulties to constrain the wavefront continuity arise when dealing with multidimensional can be biased by the choice of the path arrays. The delay within the n-D array as the delay function along the sliced path is performed by optimizing the estimate locally along the path and not globally for the overall data. The key idea underlying the least squares (LS) wavefront estimation algorithm proposed here finds its root in the analogy between the two-dimensional (2-D) wavefront estimation and the 2-D phase unwrapping. Phase unwrapping can be considered as a wavefront estimation for monochromatic signals [3]. Similar in n-D can be estito 2-D phase unwrapping, the delay mated in two steps [4], [5]. First, the differential delays between neighboring sensors are estimated by using the cross-correlation method (in 2-D phase unwrapping this is referred to as instantaneous frequency estimation), then all of these differential deby implicitly constraining lays are integrated to estimate its smoothness and continuity. The main advantage is that each step can be optimized separately for the specific application. The methods for differential delay estimation are not covered in this paper as being extensively studied in the past, the reader is referred to [6], [7] and references therein. However, if differential delays are estimated from cross correlation of signals for couples of neighboring sensors, the errors may occur in picking the lag where the cross correlation peaks. These errors, also referred to as mispicks, are caused by low SNR and limited waveform bandwidth (e.g., when dealing with narrowband waveforms, the cross correlation becomes quasi-periodic and the delay estimate can skip one or more cycles). The way to is by considering compensate the effects of these errors on the wavefront estimation from differential delay estimation as a global (weighted or unweighted) LS problem. The research has been motivated by the need to estimate the wavefront or wavefront distortion (if estimation is performed after correction with nominal shape parameters) and therefore, some of the applications are described at first (Section II). In
0196–2892/01$10.00 © 2001 IEEE
656
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 39, NO. 3, MARCH 2001
Section III, the LS algorithm for an n-D array is proposed in a common framework for regularly and irregularly spaced sensors. The LS algorithm for 2-D array is discussed in Section IV. Since the wavefront estimation depends on the differential delays, the influence of mispicks on the overall solution has been evaluated. The analysis in Section V shows that errors due to mispicks spatially decrease faster when the dimensions of the array (or equivalently the redundancy of differential delay measurements) are increased. In addition, since the mispicks are either Gaussian and correlated or non-Gaussian, the weighted LS or the iteratively reweighted LS approach are more suitable in real applications. Performances have been evaluated in Section VI by taking a reflection seismology application. repreRemark on data dimensions: Measurements sents an ( 1)-D data volume, the n-D arrangement of sources and receivers is simply indicated as the n-D array of sensors (as for the reciprocity sources and receivers can be interchanged). When sources and receivers are placed everywhere in the 6. three-dimensional (3-D) space, the array of sensors is 4 (or the Sources and receivers arranged in planes lead to data volume can be, at most, five-dimensional [5-D]), while if 2 and the data sources and receivers are along lines, then volume is 3-D. The LS method is suitable to handle n-D array of sensors provided that the signals in the neighborhood of each sensor (either source or receiver) can be accessed. In seismic 4 [8]. acquisitions, arrays are planar and II. MOTIVATIONS AND APPLICATIONS In many remote sensing applications, the properties of the propagating medium, here described by a vector that contains the model parameters, are estimated from the minimization of delays the misfit between the estimated from for locations and the modeled according traveltimes (here is the to the model parameters matrix or vector transpose) (1) The minimization strategy and the stabilization of the objective , see e.g., function mostly depends on the nonlinearity of [9]. The method proposed in this paper has been motivated by the need to reliably estimate the traveltimes when the waveform is wideband (for velocity model estimation and static corrections computation in geophysics) or narrowband (for multidimensional phase unwrapping in digital elevation mapping by SAR interferometry). The nonparametric approach has been preferred to avoid any predefined relationship between the delay estimation and the optimization (1). In this section, we revise some applications that are based on the nonparametric estimates of the delays . In other cases not discussed here, there has been an attempt to avoid the explicit delays estimation in the minimization (1). A. Statics and Traveltime Tomography In reflection seismology, the near-surface velocity anomalies are known to distort the propagating wavefield. Such distortions (statics) contribute directly to the misfocusing of deep reflec-
tors in migrated gathers. In order to correct for these effects, it is useful to limit the analysis to kinematics, thus assuming that statics can be compensated by recovering only the differential delays between neighboring traces. These time shifts can be estimated simply by cross-correlating the traces after the compensation of the nominal (or undistorted) wavefront (this processing step is usually referred to as normal moveout correction). The key to the simplicity of this approach is that all the delays are estimated independently, but this is also its main weakness, as independent estimates may lack any continuity. In order to perform (residual) static correction that minimizes the propagation of wrong delay estimates, it is necessary to exploit some of the properties of in (1). For example, the differential time de(surface consistency) in order lays can be fitted to a model to exploit data redundancy for noise attenuation [10]. Alternatively, the time shifts can be estimated by maximizing the energy of the spatial average (or stack) of the signal, thus working directly on seismic data without any model constraint [11]. Even if the latter approach leads to computationally expensive optimizations of a nonconvex objective function, it is known to cope with severe statics [12]. All of these methods exploit the redundancy implied by the vertical straight rays propagation in the near surface. This assumption is reasonable only when 1) the velocity anomaly is very shallow compared to the target-reflector depth, and 2) the velocity of shallow layers is smaller than for deeper layers. Residual static corrections may fail if the aforementioned model assumptions are violated [13]. This occurs when velocity anomalies do not give rise to a constant delay for all events in each trace and they are so severe as to induce raypath bending, or when near surface lateral velocity variations are longer than the cable and fall in the null space of the inverse problem. In all these cases, results can be improved by estimating explicitly a model for velocity anomalies by means of nonlinear traveltime [15], [16]. Furthertomography [14] or linearization of more, velocity model building cannot be avoided when the final goal is prestack depth migration of the data. Similarly to statics reflection tomography is also affected by ambiguity between velocity and depth [17]. The advantage of tomographic algorithms is that this ambiguity can be easily controlled by constraining the velocity model through any available a priori information or through any other physically meaningful constraint. In this way, the resulting model can always be made reasonable [18]. In tomographic inverse problems, model fitting (1) is carried out by comparing the estimated traveltimes with those of rays traced through the model. The traveltime estimation is usually performed by manual interpretation of prestack data, and it is largely time consuming, especially for large n-D data volume from actual acquisitions. Therefore, it would be highly desirable to have a tool for automatic traveltime estimation that is faster and possibly more accurate than “human interpretation.” B. Phase Unwrapping In synthetic aperture radar (SAR) interferometry, two (or more) backscattered fields that pertain to the same scene recorded (not necessarily at the same time) by two (or more) SAR antennas allow us to determine the position of the scattering points. The delays between different travelling paths can
BIENATI AND SPAGNOLINI: MULTIDIMENSIONAL WAVEFRONT ESTIMATION
be used to obtain the elevation map of the scene [19]. When the waveform is narrowband or monochromatic, the delays can be estimated only modulo the period of the center frequency (time dependency is (or equivalently modulo 2 ) from unnecessary and it can be dropped). In this case, the elevation map of a 2-D scene is obtained from the traveltimes only after phase unwrapping and re-scaling. 2-D phase unwrapping algorithms are based on the estimates of the differential phase (or their modifications); the phase difference for neighboring and this estimate is constrained pixels is measured modulo . The 2-D phase unwrapping to be within the interval algorithms are based on the integration of several local measurements of differential phases, each performed for the nearest neighbor of each measuring point [5]. Multiple interferograms ) phase give rise to the need of a multidimensional ( unwrapping strategy when there is the need to track the phase variations across the interferograms [20]. This problem has to deal with the simultaneous unwrapping of multiple 2-D planes, and it can be approached in a general framework as described below. Moreover, the LS approach can easily cope with irregularly sampled data. Recently, phase unwrapping on a sparse grid of stable targets has gained increasing interest (see, e.g., [21], [22]), and the use of the multidimensional approach improves considerably the reliability of the results [23]. More measurements of the same image can be carried out with the purpose of improving the SNR (i.e., the phase image shows no variation across the interferograms except for different noise realizations). In this case, the averaging of those solutions obtained by unwrapping each 2-D plane separately is equivalent to the solution that can be achieved by globally unwrapping the 3-D data volume with the constraint of having the same solution on each plane. This special case is not considered here. III. LS WAVEFRONT ESTIMATION After definition of the model (Section III-A), the LS algorithm is discussed for n-D array regardless of sensors’ positioning (Section III-B). The LS method for regularly spaced 2-D array is rederived in Section III-C by using a compact notation for structured matrixes. This is the same problem commonly solved in LS approach to 2-D phase unwrapping. Continuous-space analogy for 2-D [24] is instrumental to derive the properties of LS method for n-D array (Section IV). A. Data model 1)-D data volume Let the ( layed waveform
be described as a de-
(2) has been (approxiwhere the target event is any mately) windowed from the whole data volume, spatially and temporally structured noise (i.e., all the interfering influences only the events), and the amplitude variations for the target reflection (or for any SNR. The waveform compound reflection) need not be known and can also be slowly , with changing in space: . In other words, for any neighboring
657
point , the normalized correlation (or the coherence between the two waveforms)
(3)
1. The delay function is should be assumed to be a continuous and regular surface (up to the second order derivatives), and the continuity also implies that has to be single valued (without any triplications or faults). In the presence of diffractions, it may be necessary to backpropagate the wavefield (migration), estimate the delay , or in the migrated domain, and then demior , see e.g., [25]. grate the estimated delay B. LS algorithm Let arg
(4)
and estibe the differential delay between the signals in mated by using the cross-correlation method, the unknown travcan be estimated by solving, in the LS sense, eltime surface for the equality (5) measurements. The estimate within a predefinite subset of is thus given by the solution of the linear system (6) equations for unknowns (usually ), is the of delays, and is a vector consisting of the difvector of ferential delays between all the couples of traces useful for the differential delay estimation. is a sparse matrix, each row has only two nonzero elements that corresponds to the differential delay between the two sensors. The LS solution of the system (6) is known to be given by (7) and it can be computed efficiently by means of iterative methods that exploit the sparsity of (e.g., successive over-relaxation, preconditioned conjugate gradient, or QR decomposition; see [26] or [27] for a detailed discussion). The system (6) and its LS solution (7) represent the general formulation of the problem for n-D arrays, as the regularity of the array geometry is not required. Only the selection of the couples of sensors for the set of differential delay estimates is needed. The preferred choice is for those couple of traces where the waveform is slowly changing and the signal is mostly spatially-coherent, i.e., 1. In practice, the subset is lower than the sensors: ( 1)/2. overall number of couples for is in order. The A discussion on the structure of matrix can be carried out by sensors’ pairing to obtain the matrix considering a neighboring selection criteria. The subset of differential delays can be estimated by looking for the sensors pairs
658
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 39, NO. 3, MARCH 2001
and so that the geometrical distance , for . Depending on the array geometry, the number some value of measurements is linearly increasing with the number of sen, where is the number of neighbors to each sors: sensor. For regularly spaced sensors, the selection criteria for (the system surrounding sensors is trivial and it yields 2). In addition, when the ratio (6) is overdetermined only if between the number of constraints and the unknowns is constant within a large support (or in the whole support as for regularly spaced array), the solution is expected to have predictable properties. For an irregular acquisition geometry, the problem of determining the nearest neighbor of each source and/or receiver can be solved by using any of the known methods. However, here we preferred the use of the Delaunay triangulation algorithm [28] even if the sensors’ pairing algorithm is not critical as far as the selection criteria is reasonable. The structure of matrix depends only on the sensors ordering. The properties and the limitations of the LS approach can be evaluated only for some simple (but useful) cases with regular array geometry: the 2-D array and the n-D array. Most of the analytic properties of the LS approach can be more easily derived when considering the continuous-space wavefront (not sampled). This will be discussed in Section V.
is a small constant value that stabilizes the iterations. The . LS solution (7) represents the first iteration, i.e., It can be shown that this choice of the weighting matrix leads the IRLS to have performances similar to the L1-norm [29] and thus, IRLS turns out to be quite effective in reducing the influence of large mispicks in differential delays . IV. WAVEFRONT ESTIMATION FOR 2-D ARRAY For a uniformly spaced 2-D array, the differential delay for neighboring traces can be evaluated along the two axes here redenote the maferred to as and . Let R and trix of (unknown) delays, R the matrices of differential dethe lays for neighboring traces (or gradients), and shift matrix, i.e., where is the Kronecker 1 for 0, and 0 for 0). The matrix operdelta ( ators to compute the gradients from the matrix can be written and , where as is the unitary matrix. Let and be the differential delay matrix augmented by or zeros, respectively. The delays can be oba vector of tained from the following linear relationships on the gradients: (9)
C. Weighted and iteratively re-weighted LS The equality (5) can be weighted according to any predefinite criteria. For a positive definite weighting matrix (not necessarily diagonal), the minimization of the weighted leads to the following solution: norm . In addition to the weighting matrix , it is possible, when the Gaussian model holds true, to introduce the a priori information through the covariance matrix of traveltimes. Since these matrices are often unknown, feasible approximations are sometimes needed. The covariance matrix of can be chosen so that the dispersion of each measurement is proportional to the geometrical distance , or equivis diagonal with entries set as . This was alently adopted in the seismic applications discussed in Section II-A. with 0 1 can be used to Diagonal matrix pair the sensors according to any selection criteria based on the stationarity of either the wavefield or the SNR. In principle, at the beginning, all possible sensor pairs can be taken into ac( 1)/2. Then, after some iterations, terms are count and 0. In this latter case, downweighted or can be set to the corresponding equations are simply removed from system (6). can also be computed from The weighting matrix , and the LS solution (7) can be the residuals solved iteratively. This technique is also known as the iteratively be the estireweighted least squares (IRLS). In practice, let mated wavefront at the th iteration according to solution of the , weighted LS system for the next iteration depends on the the weight matrix , and it is diagonal overall residual
(8)
vec According to the equality vec for the vec( ) operator that stacks all columns of a matrix into one single column vector and some of the properties of Kronecker products [30], the system of (9) can be vec , rewritten similar to system (6) by defining vec vec and . , and the equaThe number of unknowns is 2 . For tions are matrix regularly sampled 2-D array, the is structured and sparse similar to Hunt [31]. The normal equations can be rewritten as (10) From the property of the Kronecker products, the normal equations [(10)] are equivalent to the discrete version of the Poisson’s equation comprehensive of the Neumann boundary (or ) conditions [3]. For instance, the matrix is the discrete version of the second order derivative along (or ). is identically zero, as it is associOne eigenvalue of ated to the eigenvector corresponding to constant shifts (this is expected as the LS wavefront estimation is based on the differential delays). For 2-D array, the error analysis becomes easy if the estimation problem is reconsidered in continuous space. Assume that we have estimated the gradient field (where and are the unit norm vectors for the and axes, respectively) of the delay function within the region of support . The wavefront can be estimated from the and , the minimization of the misfit between
BIENATI AND SPAGNOLINI: MULTIDIMENSIONAL WAVEFRONT ESTIMATION
gradient of (here ). The solution of the variational problem is provided by the integration of the Poisson’s equation (11) and are the Laplace and the divergence opwhere erators, respectively. The numerical solution of (11) with Neumann boundary conditions obtained by approximating partial derivatives with finite differences leads to (10). When (11) is used for the wavefront estimation, it has a nonempty data null space. This null space contains those components of the differential delays that cannot be predicted by the model, and it represents a residual misfit after the minimum has been achieved. This is due to the choice of describing the gradients field by using the scalar function only. According to the Helmoltz’s theorem, any vector can be represented by the superposition of a scalar field and a vector potential potential (12) cannot give rise to where the continuous surface itself , which is accounted for by the rotational component of the vector potential. This component must necessarily be part of 0 for any . the data null space as The rotational component is due to those values in the estimated gradients that conflict with the continuity of the solution caused either by mispicks and/or by not complying with the continuity hypothesis (e.g., when the surface is faulted or represents a linear transformamultivalued). Since matrix tion, the rotational component is also given by the projection of the measurements onto the null space of (see [3] for further details). When dealing with irregularly spaced array, the evaluation of the projection matrix is the only way to establish the sources of error. In any case, similar to 2-D phase unwrapping, the rotational component field represents a useful indicator of the quality of the cross-correlation maxima picking and of the match of the delay surface with the continuity hypothesis. V. ADVANTAGES OF WAVEFRONT ESTIMATION MULTIDIMENSIONAL ARRAYS
IN
In principle, n-D array can be sliced to extract several 2-D arrays and have simpler wavefront estimation algorithms. Here we demonstrate, by simple arguments, that this approach is far from being an optimum one except for trivial cases. A multidimensional approach for n-D array is expected to have a better constraint as, in the simpler case of regularly sampled arrays, the equations (one delays are obtained by solving a system of unknowns. Therefore, equation for each dimension) for the 3), a higher redundancy of mulfor ( 1)-D data volumes ( tidimensional algorithms guarantees better performances when compared to iterated estimation on lower dimensional subvolumes obtained by slicing the global volume. The advantages for n-D array can be evaluated in terms of error propagation by considering the continuous space case. Let contain only the estimated gradients one single incorrect measure due to a mispick (for instance in along one component only, say :
659
), the wavefront error caused by follows from the integration of the partial differential equation (13) , it is derived from (11) as for the true wavefront . On each dimension, the finite differences approximation splits the mispick into two neighboring terms of of opposite sign: a dipole. Since the driving term in (13) is a dipole, the solution is the combination of the responses of two pulses having opposite signs that decays faster than each pulse response. can be evaluated first by evaluating the The decaying of to one pulse: . response of for According to Gauss’s law, the amplitude of decays as the inverse of the surface of the [32]. For a dipole having unit width n-D sphere: spacing the decays is the combination of both terms (when is ) larger than the dipole width, in this case where even (14) odd or equivalently [by assuming that for radial]
the field (14) is
(15) In conclusion, the advantage of the multidimensional traveltime 2 is the lower spreading of those errors estimation for that lead to conflicting values in the gradient (also referred to as rotational sources in 2-D phase unwrapping [4]), these errors are mainly due to mispicks. This qualitative proof that mixes continuous and discrete doin (15) is given mains is far from being a rigorous one as consists of by the solution of (13) when the forcing term one pulse and it is measured far from the pulse-location. For the discrete-space the error behavior depends on the numerical solution of (7) that cannot be recursively computed as its support is not wedge shaped, the results are shown in Fig. 1 for 2-D and 3-D arrays. In both cases the mispick occurs on , the 3-D response is shown only on one plane as it is invariant by rotation. According to the relationship (15) the error (asymptotically) defor 3-D arrays. Moreover, also the cays as 1 for 2-D and 1 peak amplitude decreases as the number of spatial dimensions increases (from 1/4 for 2-D to 1/6 for 3-D arrays). When the errors are far from being Gaussian distributed (as for errors that arise from mispicks), the approximate L-1 norm approach or IRLS is advisable. VI. EXAMPLES AND PERFORMANCES FROM REFLECTION SEISMOLOGY In reflection seismology, the space-axes represent any suitable combination of space coordinates from the acquisition ge-
660
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 39, NO. 3, MARCH 2001
Fig. 1. Mispick response (x) for wavefront estimation in two-dimensional (2-D) array (left: the response has been magnified by a factor of 4) and in three-dimensional (3-D) array (right: the response has been magnified by a factor of 6).
Fig. 2.
Source-receiver geometry for 3-D land acquisiton (SEG/EAEG 3-D Overthrust B-3 acquisition subset).
ometry (source and receiver positions, midpoint and offset, in 2-D experiments; receiver positions in a common shot gather, offset-azimuth or any other suitable combination for 3-D acquisitions). Even if the LS wavefront estimation algorithm is general, here we tailor the method for 3-D and four-dimensional (4-D) data volumes. A short description of the 3-D acquisition geometry will be provided at first (for a detailed descriptions of source/receiver arrangements the reader is referred to [8]). The way how the source/receiver arrangements sample the wavefield depends on the specific acquisition templates [33]. denotes the coordinate associated to the In this section ) those associated to the receivers. For field sources and ( acquisitions the five-dimensional (5-D) data volume consists of sources and receivers, the overall 4-D array samples the wavefield into locations. The 3-D land acquisition geometry considered here has multiple cross-spreads (each cross-spread is a 3-D data volume as sources and receivers moves along orthogonal lines), as shown , in Fig. 2. The 4-D volume of the gathered data becomes and . The data considered for here
the examples proposed here are from the SEG/EAEG 3-D Overthrust B-3 acquisition subset (dataset available at [34], see [35] for a description of the experiment). According to Fig. 2, we selected eight consecutive source-receiver lines (cross-spreads) that correspond either to eight traveltime surfaces each obtained for a given or to an overall traveltime volume . The main feature of this acquisition is the nonuni325 m 50 m form spacing of receiver lines where for , 8. The sensors are uniformly spaced along the receiver lines (spacing is 50 m). The reflections of the target event arise from a nonflat reflector embedded in complex interfering backscattered wavefields (data and the estimated volume are in Fig. 3. This subset is realistic enough to be used here to compare the performances of the LS method for varying SNR, increasing dimensions, and sparse/irregular acquisitions. IRLS is also evaluated. The advantages of the wavefront estimation for multidimensional arrays is shown in Fig. 4, which quantifies, for the data 2 receiver lines in 350 m volume obtained from 400 m, the better immunity to noise when dealing with and
BIENATI AND SPAGNOLINI: MULTIDIMENSIONAL WAVEFRONT ESTIMATION
661
Fig. 3. Slices from the subset of SEG/EAGE three-dimensional (3-D) Overthrust B-3 data. 3-D data volume d(x ; x = 325m; y ; t) sliced along three axes [planes in gray-scale: d(x ; x = 325 m, y ; t = 0.85 s), d(x = 2000 m, x = 325 m, y ; t), and d(x ; x = 325 m, y = 900 m, t); white lines: signals in four receiver-locations] and the estimated delays (x ; x ; y ) for the first cross-spread (x = 325 m). Wavefront correction has been applied for visual display to flatten the nearly hyperbolic wavefront.
Fig. 4. Root mean square error (RMSE) of estimated delays versus noise percentage for 2-D array (or 3-D data volume) and 3-D array (or 4-D data volume). Experimental results after 50 runs (N = 2, N = 120, N = 31, N = 1).
4-D data volumes compared to 3-D. The comparison is given in terms of root mean squared error (RMSE) of estimated traveltime for increasing levels of (white Gaussian) noise superimposed to the data set. The four surfaces (estimated separately on each 3-D data volume) and the volume (estimated globally from the 4-D data volume) obtained in the absence of Gaussian noise have been used as references to evaluate the RMSE after the noise has been added. The standard deviation of the noise level was normalized to the root mean squared value of the signal as this can be more easily related to the way the noise level is
Fig. 5. Wavefront estimation (x ; x ; y ) for 3-D array (cross-spread acquisition). The 3-D delay surface is estimated by using the full 4-D data volume are plotted for varying x = {325 m, 625 m, 925 m, 1225 m} (view rotated by 90 with respect to Fig. 3).
inferred from experimental measurements. In Fig. 6 the same comparison between the strategies that deal with 3-D and 4-D data volumes is shown for 95% and 160% of noise. In any case, for noise below 120% the RMSE is approx. below the sampling interval (8 ms). The advantages of higher dimensions, mainly for higher noise level, are justified by the larger redundancy. In addition, the advantages of the IRLS approach for large noise is a consequence of the interaction between the noise and the cycle-skips for low SNR.
662
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 39, NO. 3, MARCH 2001
=
Fig. 6. Slices of data volume (same as Fig. 3) for the first cross-spread (x 325 m) and estimated wavefront for increasing additive noise level. Left column: the wavefront is estimated in 2-D array. Right column: the wavefront is estimated in 3-D array (the same slice as for 2-D array is shown in right column).
When the estimation is carried out through subvolumes slicing, the estimated wavefronts can be mistied as each estimate cannot be easily tied together with the others (in some cases, the slicing strategy can be arranged to have several common points among different slices). For example, in acquisition comprising multiple cross-spreads (see Fig. 2) each traveltime surface can be estimated from each cross-spread can be (2-D array) or the overall delay volume estimated from the whole data set. Once estimated, the volume can be sliced into surfaces if necessary, but each of these surfaces are naturally tied by multidimensional estimation and they are not floating as for traveltime estimation after slicing (where the relative shift of each solution is arbitrary). This is another reason to exploit the full dimensionality of the data set. Fig. 5 shows the four traveltime surfaces (for 325 m, 625 m, 925 m, 1225 m) obtained by the 4-D algorithm from the mentioned subset of SEG/EAEG 3-D Overthrust. Due to the slope of the target reflector the surfaces have a different shape but 4-D estimation has correctly tied the delay surfaces shows an increasing where the reflector is flat and delay when the receiver lines move up-dip. The same data set considered in Fig. 6 has been randomly decimated 2:3 and 2:5 in order to simulate an acquisition with sparse source-receivers randomly located. The delay function
has been estimated on a single decimated cross-spread (2-D array). Fig. 7 shows that even in absence of additive noise the reduction of the number of traces increases the error in the estimated wavefront due to a lower coherence (3) between the waveform across the array caused by a larger distance between neighboring sensors. Obviously the results worsen when white Gaussian noise is added to the data (Fig. 7, right). In order to keep the same redundancy for different decimations, the differential delays have been estimated on a number of pairs so that the ratio is approximately con(or equations) stant. The advantage of the IRLS method is shown in Fig. 8 for 4-D data volumes and increasing noise level. Here the first 3 receiver lines have been considered having sensors’ {325 m, 375 m, 625 m}, respectively. The spacing of greater distance between the second and the third line (250 m) implies that the slowly changing waveform hypothesis is not fully satisfied, and it increases the probability of cycle-skip occurrence in cross-correlation maxima picking. Even if the large array spacing makes the advantages of the multidimensional approach rather negligible, the wavefronts are tied consistently (e.g., Fig. 5). In addition, the interaction between high noise level and large errors (mainly due to cycle-skips in differential delay estimation) makes the IRLS method an attractive tool to handle low SNR data.
BIENATI AND SPAGNOLINI: MULTIDIMENSIONAL WAVEFRONT ESTIMATION
663
=
Fig. 7. Wavefront estimation in 2-D array for decimated data extracted from the first cross-spread (x 325 m). Left: wavefront estimation increasing decimation rate (2:3 and 2:5, noise free data). Right: wavefront estimation for noisy data (95% of Gaussian noise) and increasing decimation rate.
Fig. 8. RMSE of estimated delays for 3-D array (irregular spacing for x ) versus noise percentage for LS approach and IRLS. Experimental results after 50 runs ( 3, 120, 31, 1).
N = N =
N =
N =
VII. CONCLUSIONS This paper has shown that a least squares approach (or any modifications such as weighted LS or iteratively reweighted
LS) can be used for the nonparametric estimate of a delay (hyper)volume from the differential delays for neighboring sensors. Cross-correlation maxima picking is used as differential delays estimators. The errors due to peak selection ambiguity of cross-correlator estimator (mispick) have been evaluated and their propagation has been proved to be limited by the LS method that exploits the spatial continuity of the wavefront over the multidimensional space. The robustness against noise and mispicks improves with the dimensionality of the data set. The L1-norm approach is advisable when the mispicks are more likely (e.g., when the sensor’s spacing is too large and the waveform changes across the array). In the proposed method, the constraint that the differential estimates are only for neighboring sensors has been relaxed. Therefore, it can be extended to the differential delay estimates for any sensors’ configuration (and thus for irregularly spaced sensors) and any waveform (wideband or pulsed waveforms and narrowband or phase unwrapping problems). ACKNOWLEDGMENT The authors would like to thank F.Rocca for some stimulating discussions on the application of the multidimensional approach to geophysics.
664
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 39, NO. 3, MARCH 2001
REFERENCES [1] U. Spagnolini and S. Grion, “Shape parameter estimation of wavefronts with known waveform,” Signal Process., vol. 68, pp. 233–257, Aug. 1998. [2] D. H. Johnson and D. E. Dudgeon, Array Signal Processing: Concepts and Techniques. Englewood Cliffs, NJ: Prentice-Hall, 1993. [3] U. Spagnolini, “Non-parametric narrow-band wavefront estimation from wavefront gradients,” IEEE Trans. Signal Processing, vol. 47, pp. 3116–3121, Nov. 1999. , “2-D phase unwrapping and instantaneous frequency estimation,” [4] IEEE Trans. Geosci. Remote Sensing, vol. GRS-33, pp. 579–589, May 1995. [5] D. C. Ghiglia and M. D. Pritt, Two-Dimensional Phase Unwrapping. New York: Wiley, 1998. [6] IEEE Trans. Acoust., Speech, Signal Processing, vol. ASSP-29, Pt. II, June 1981. [7] G. C. Carter, “Coherence and time delay estimation,” Proc. IEEE, vol. 75, pp. 236–255, Feb. 1987. [8] D. G. Stone, Designing Seismic Surveys in Two and Three Dimensions. Tulsa, OK: Soc. Expl. Geophys., 1994. [9] J. A. Scales, P. Docherty, and A. Gersztenkorn, “Regularization of nonlinear inverse problem: Imaging the near-surface weathering layer,” Inv. Probl., vol. 6, pp. 115–131, 1990. [10] R. A. Wiggins, K. L. Larner, and +R. D. Wisecup, “Residual statics analysis as a general linear inverse problem,” Geophysics, vol. 41, no. 5, pp. 922–938, Oct. 1976. [11] J. Ronen and J. F. Claerbout, “Surface consistent residual statics estimation by stack-power maximization,” Geophysics, vol. 50, no. 12, pp. 2759–2767, Dec. 1985. [12] D. H. Rothman, “Nonlinear inversion, statistical mechanics, and residual statics estimation,” Geophysics, vol. 50, no. 12, pp. 2784–2796, Dec. 1985. [13] D. Palmer, Refraction Seismics. London, U.K.: Geophysical Press, 1986. [14] T. N. Bishop, K. P. Bube, R. T. Cutler, R. T. Langan, P. L. Love, J. R. Resnick, R. T. Shuey, D. A. Spindler, and H. W. Wyld, “Tomographic determination of velocity and depth in laterally varying media,” Geophysics, vol. 50, pp. 903–923, June 1985. [15] E. Loinger, “A linear model for velocity anomalies,” Geophys. Prospect., vol. 31, pp. 98–118, Feb. 1983. [16] F. Rocca and U. Spagnolini, “A simple interpretation of reflection tomography,” in Proc. Eur. Assoc. Exploration Geophysicists (EAEG) Meeting, vol. 1, Glasgow, U.K., 1995, pp. D2–D3. [17] C. Stork, “Singular value decomposition of the velocity-reflector depth tradeoff, Part 2: High-resolution analysis of a generic model,” Geophysics, vol. 57, pp. 933–943, July 1992. [18] C. Stork and R. W. Clayton, “Using constraints to address the instabilities of prestack velocity analysis,” Geophysics, vol. 57, pp. 404–419, Mar. 1992. [19] H. A. Zebker and R. M. Goldstein, “Topographic mapping from interferometric synthetic aperture radar observations,” J. Geophys. Res., vol. 91, pp. 4993–4999, 1986. [20] A. Ferretti, C. Prati, F. Rocca, and A. Monti Guarnieri, “Multi-baseline SAR interferometry for automatic DEM reconstruction,” in Proc. 3rd ERS Symp., Florence, Italy, 1997. [21] M. Costantini and P. Rosen, “A generalized phase unwrapping approach for sparse data,” in Proc. Int. Geosci. Remote Sensing Symp., Hamburg, Germany, 1999, pp. 267–269. [22] M. Eineder and J. Holzner, “Phase unwrapping of low coherence differential interferograms,” in Proc. Int. Geosci. Remote Sensing Symp., Hamburg, Germany, 1999, pp. 1727–1730.
[23] A. Ferretti, C. Prati, and F. Rocca, “Non-linear subsidence rate estimation using permanent scatterers in differential SAR interferometry,” IEEE Trans. Geosci. Remote Sensing, vol. 38, pp. 2202–2212, Sept. 2000. [24] H. Takajo and T. Takahashi, “Noniterative method for obtaining the exact solution for the normal equation in least-squares phase estimation from the phase difference,” J. Opt. Soc. Amer. A, vol. 5, pp. 1818–1827, 1988. [25] J. A. C. Jacobs, F. Delprat-Jannaud, A. Ehinger, and P. Lailly, “Sequential migration-aided reflection tomography: A tool for imaging complex structures,” in 62nd Annu. Int. Meeting SEG, Expanded Abstract, New Orleans, LA, 1992. [26] G. H. Golub and C. F. Van Loan, Matrix Computations, 2nd ed. London, U.K.: Johns Hopkins Univ. Press, 1989. [27] G. Meurant, Computer Solution of Large Linear Systems. Amsterdam, The Netherlands: Elsevier, 1999. [28] F. Preparata and M. Shamos, Computational Geometry—An Introduction. Berlin, Germany: Springer-Verlag, 1985. [29] J. A. Scales and A. Gersztenkorn, “Robust method in inverse theory,” Inv. Probl., vol. 4, pp. 1071–1091, 1988. [30] A. Graham, Kronecker Products and Matrix Calculus. New York: Wiley, 1981. [31] B. R. Hunt, “Matrix formulation of the reconstruction of phase values from phase differences,” J. Opt. Soc. Amer. A, vol. 69, no. 3, pp. 393–399, Mar. 1979. [32] R. Courant and T. Hilbert, Methods of Mathematical Physics. New York: Wiley. [33] G. J. O. Vermeer, Seismic Wavefield Sampling. Tulsa, OK: Soc. Expl. Geophys., 1990. [34] SEG/EAEG 3-D Modeling data access site: http://archiv.llnl.gov/SSD, . [35] F. Aminzadeh, N. Burkhard, P. Lailly, F. Rocca, and K. Wyatt, “Progress Rep. from the SEG/EAEG 3-D modeling commitee,”, vol. 13, 1994.
Nicola Bienati received the Dott. Ing. degree (cum laude) in telecommunications and the Ph.D. degree in electronics and telecommunications from the Politecnico di Milano, Milano, Italy, in 1995 and 2000, respectively. His main research interests are seismic reflection tomography, horizon picking in multidimensional data volumes, blind deconvolution, and blind sources separation.
Umberto Spagnolini (M’99) received the Dott. Ing. Elett. degree in telecommunications (cum laude) from the Politenico di Milano, Milano, Italy, in 1988. Since 1988, he has been with the Dipartimento di Elettronica e Informazione, Politenico di Milano, where he has held the position of Associate Professor of Digital Signal Processing since 1998. His primary research (and applications) focuses on array processing and wavefield interpolation (mobile communication and geophysics), inverse problems (ground penetrating radar) and parameter estimation (2-D phase unwrapping for SAR), and non-Gaussian EMI reduction. Dr. Spagnolini is a member of the SEG and EAGE and serves as an Associate Editor for IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING. He was awarded the Associazione Elettrotecnica ed Elettronica Italiana (AEI) Award and the Van Weelden Award of EAGE, both in 1991, and received the Best Paper Award from EAGE in 1998.
