New advances in three-dimensional controlled-source electromagnetic inversion

New techniques for improving both the computational and imaging performance of the three-dimensional (3-D) electromagnetic inverse problem are presented. A non-linear conjugate gradient algorithm is the framework of the inversion scheme. Full wave equation modelling for controlled sources is utilized for data simulation along with an efficient gradient computation approach for the model update. Improving the modelling efficiency of the 3-D finite difference (FD) method involves the separation of the potentially large modelling mesh, defining the set of model parameters, from the computational FD meshes used for field simulation. Grid spacings and thus overall grid sizes can be reduced and optimized according to source frequencies and source–receiver offsets of a given input data set. Further computational efficiency is obtained by combining different levels of parallelization. While the parallel scheme allows for an arbitrarily large number of parallel tasks, the relative amount of message passing is kept constant. Image enhancement is achieved by model parameter transformation functions, which enforce bounded conductivity parameters and thus prevent parameter overshoots. Further, a remedy for treating distorted data within the inversion process is presented. Data distortions simulated here include positioning errors and a highly conductive overburden, hiding the desired target signal. The methods are demonstrated using both synthetic and field data.     

Non-linear inversion using general measures of data misfit and model structure

We investigate the use of general, non- l ₂ measures of data misfit and model structure in the solution of the non-linear inverse problem. Of particular interest are robust measures of data misfit, and measures of model structure which enable piecewise-constant models to be constructed. General measures can be incorporated into traditional linearized, iterative solutions to the non-linear problem through the use of an iteratively reweighted least-squares (IRLS) algorithm. We show how such an algorithm can be used to solve the linear inverse problem when general measures of misfit and structure are considered. The magnetic stripe example of Parker (1994 ) is used as an illustration. This example also emphasizes the benefits of using a robust measure of misfit when outliers are present in the data. We then show how the IRLS algorithm can be used within a linearized, iterative solution to the non-linear problem. The relevant procedure contains two iterative loops which can be combined in a number of ways. We present two possibilities. The first involves a line search to determine the most appropriate value of the trade-off parameter and the complete solution, via the IRLS algorithm, of the linearized inverse problem for each value of the trade-off parameter. In the second approach, a schedule of prescribed values for the trade-off parameter is used and the iterations required by the IRLS algorithm are combined with those for the linearized, iterative inversion procedure. These two variations are then applied to the 1-D inversion of both synthetic and field time-domain electromagnetic data.     

Least-squares fitting of marine seismic refraction data

Summary. An iterative procedure is presented for fitting waveform data from a marine seismic refraction experiment. The wavefunction from the explosive source is known and the crustal structure is refined using the damped least squares procedure. The damping parameter serves the dual purpose of stabilizing an under-constrained inversion and improving the linearity by suppressing high frequencies. The synthetic seismograms and their model differentials are calculated using the WKBJ seismogram method. Both the synthetic seismograms and the linear algebra are sufficiently straightforward that the computations can be performed on an array processor. The inversion procedure is then sufficiently rapid that interactive computations are possible. The technique is illustrated using the FF2 refraction data from the Fanfare cruise of the Scripps Institution of Oceanography. These data had been interpreted previously by trial-and-error using the reflectivity method. Starting from two different, simple models, the inversion procedure obtains essentially one unique model. Its features are very similar to the previous model.     

Simulated envelopes of non-isotropically scattered body waves as compared to observed ones: another manifestation of fractal heterogeneity

Envelopes of scalar waves are simulated at various distances from an instant point source embedded in a random uniformly scattering medium by means of direct Monte-Carlo modelling of wave-energy transport. Three types of scattering radiation pattern ('indicatrix') are studied, for media specified by (1) a Gaussian autocorrelation function of inhomogeneities, (2) a power-law ('fractal', k -^α) inhomogeneity spectrum and (3) the mix of case (1) and the isotropic indicatrix (very small + large inhomogeneities). We look for a model that can qualitatively reproduce the two most characteristic features of real S-wave envelopes of near earthquakes, namely (1) the broadening of the 'direct' wave group with distance and (2) the monotonously decaying shape of the coda envelope that does not deviate strongly from that expected in the isotropic scattering case. Both properties are observed for any band over a wide frequency range (1-40 Hz). The well-studied isotropic scattering model realistically predicts the appearance of codas but fails to predict pulse broadening. The model of large-scale inhomogeneity realistically predicts the mode of pulse broadening but fails to predict codas. We have found that, for a particular frequency band, within each class of inhomogeneity studied, both requirements can be qualitatively satisfied by a certain choice of parameters. In the Gaussian-ACF case, however, this match can be obtained only for a narrow frequency range. In contrast, the fractal case (with a value of exponent a of about 3.5-4) reproduces qualitatively the observed wide-band behaviour, and we consider it a reasonable representation of the gross properties of the earth medium.     

Three-dimensional lithospheric modelling at NORSAR: linearity of the method and amplitude variations from the anomalies

Summary. A modification of the Aid et al . technique for three-dimensional lithospheric modelling is used to find smoothly varying models for the P -wave velocity structure beneath NORSAR. The method includes ray tracing and calculation of geometrical spreading in the anomalies. The results of linear inversion of the travel-time data compare well with those of previous investigators. The assumption of linearity, which removes the need to ray trace through the anomalies, is tested with iterative solutions for both synthetic and real data. A model with an rms velocity perturbation of 3 per cent, extending to 120 km depth, is found to be reasonably linear. In fact the procedure leads to two models which satisfy the same amount of the real data but which differ by far more than the standard errors. However, these differences are not significant once the imperfect resolution is accounted for by using the total estimation error of the stochastic inverse.
The depth of major anomalies appears to be greater than the array diameter and is therefore not well constrained. Comparing the geometrical spreading produced by these models with the amplitude variations observed at the array indicates that structure deeper than 120 km but shallower than 200 km makes an important contribution to the observations. None of the models used can produce variations as large as those in the amplitude data. For deep, essentially two-dimensional, anomalies the fit to these data is much better for sources to the NE of the array than for sources in other quadrants.     

Sequential integrated inversion of refraction and wide-angle reflection traveltimes and gravity data for two-dimensional velocity structures

A new algorithm is presented for the integrated 2-D inversion of seismic traveltime and gravity data. The algorithm adopts the 'maximum likelihood' regularization scheme. We construct a 'probability density function' which includes three kinds of information: information derived from gravity measurements; information derived from the seismic traveltime inversion procedure applied to the model; and information on the physical correlation among the density and the velocity parameters. We assume a linear relation between density and velocity, which can be node-dependent; that is, we can choose different relationships for different parts of the velocity–density grid. In addition, our procedure allows us to consider a covariance matrix related to the error propagation in linking density to velocity. We use seismic data to estimate starting velocity values and the position of boundary nodes. Subsequently, the sequential integrated inversion (SII) optimizes the layer velocities and densities for our models. The procedure is applicable, as an additional step, to any type of seismic tomographic inversion.
We illustrate the method by comparing the velocity models recovered from a standard seismic traveltime inversion with those retrieved using our algorithm. The inversion of synthetic data calculated for a 2-D isotropic, laterally inhomogeneous model shows the stability and accuracy of this procedure, demonstrates the improvements to the recovery of true velocity anomalies, and proves that this technique can efficiently overcome some of the limitations of both gravity and seismic traveltime inversions, when they are used independently.
An interpretation of field data from the 1994 Vesuvius test experiment is also presented. At depths down to 4.5 km, the model retrieved after a SII shows a more detailed structure than the model obtained from an interpretation of seismic traveltime only, and yields additional information for a further study of the area.     

青藏高原地区的大气浑浊度特征

本文通过比较Angstrom大气浑浊度系数的各种计算方法,得出了适合青藏高原地区的大气浑浊度系数的计算方法。据此利用1982年8月至1983年7月青藏高原地区部分观测站的太阳直接辐射资料计算了大气浑浊度系数,分析了它的变化特征。同时利用常规气象要素拟合得出了计算大气浑浊度系数的经验公式,计算和分析了青藏高原地区大气浑浊度系数的时间和空间分布。     

基于最大似然估计的海面风场反演算法研究

最大似然估计法被认为是海面风场反演的最佳方法,目前被用来处理SeaWinds散射计数据。风矢量求解算法是风场反演算法的核心内容。最大似然法的目标函数形式决定了在风场反演过程中必须采用数值方法求解风矢量,而传统数值求解方法运算量大。该文对最大似然估计的风场反演方法的基本原理和具体过程进行探讨,根据其目标函数的一般分布特征,提出一种较为高效的数值风矢量搜索算法。用SeaWinds散射计的12A实测数据和相应的L2B数据验证了该算法的可行性。     

A search for the ratio between gravity variation and vertical displacement due to a surface load

A data space approach to magnetotelluric (MT) inversion reduces the size of the system of equations that must be solved from M × M , as required for a model space approach, to only N × N , where M is the number of model parameter and N is the number of data. This reduction makes 3-D MT inversion on a personal computer possible for modest values of M and N . However, the need to store the N × M sensitivity matrix J remains a serious limitation. Here, we consider application of conjugate gradient (CG) methods to solve the system of data space Gauss–Newton equations. With this approach J is not explicitly formed and stored, but instead the product of J with an arbitrary vector is computed by solving one forward problem. As a test of this data space conjugate gradient (DCG) algorithm, we consider the 2-D MT inverse problem. Computational efficiency is assessed and compared to the data space Occam's (DASOCC) inversion by counting the number of forward modelling calls. Experiments with synthetic data show that although DCG requires significantly less memory, it generally requires more forward problem solutions than a scheme such as DASOCC, which is based on a full computation of J .     

Inversion of seismic attributes for velocity and attenuation structure

We have developed an inversion formuialion for velocity and attenuation structure using seismic attributes, including envelope amplitude, instantaneous frequency and arrival times of selected seismic phases. We refer to this approach as AFT inversion for amplitude, (instantaneous) frequency and time. Complex trace analysis is used to extract the different seismic attributes. The instantaneous frequency data are converted to t * using a matching procedure that approximately removes the effects of the source spectra. To invert for structure, ray-perturbation methods are used to compute the sensitivity of the seismic attributes to variations in the model. An iterative inversion procedure is then performed from smooth to less smooth models that progressively incorporates the shorter-wavelength components of the model. To illustrate the method, seismic attributes are extracted from seismic-refraction data of the Ouachita PASSCAL experiment and used to invert for shallow crustal velocity and attenuation structure. Although amplitude data are sensitive to model roughness, the inverted velocity and attenuation models were required by the data to maintain a relatively smooth character. The amplitude and t * data were needed, along with the traveltimes, at each step of the inversion in order to fit all the seismic attributes at the final iteration.     

Smooth inversion for ground surface temperature histories: estimating the optimum regularization parameter by generalized cross-validation

The inversion of recent borehole temperatures has proved to be a successful tool to determine ancient ground surface temperature histories. To take into account heterogeneity of thermal properties and their non-linear dependence on temperature itself, a versatile 1-D inversion technique based on a finite-difference approach has been developed. Regularization of the generally ill-posed problem is obtained by an appropriate version of Tikhonov regularization of variable order. In this approach, a regularization parameter has to be determined, representing a trade-off between data fit and model smoothness. We propose to select this parameter by generalized cross-validation. The resulting technique is employed in case studies from the Kola ultradeep drilling site, and another borehole from northeastern Poland. Comparing the results from both sites corroborates the hypothesis that subglacial ground surface temperatures as met in Kola often are much higher than the ones in areas exposed to atmospheric conditions (Poland).     

Three-dimensional magnetotelluric inversion using non-linear conjugate gradients

We have formulated a 3-D inverse solution for the magnetotelluric (MT) problem using the non-linear conjugate gradient method. Finite difference methods are used to compute predicted data efficiently and objective functional gradients. Only six forward modelling applications per frequency are typically required to produce the model update at each iteration. This efficiency is achieved by incorporating a simple line search procedure that calls for a sufficient reduction in the objective functional, instead of an exact determination of its minimum along a given descent direction. Additional efficiencies in the scheme are sought by incorporating preconditioning to accelerate solution convergence. Even with these efficiencies, the solution's realism and complexity are still limited by the speed and memory of serial processors. To overcome this barrier, the scheme has been implemented on a parallel computing platform where tens to thousands of processors operate on the problem simultaneously. The inversion scheme is tested by inverting data produced with a forward modelling code algorithmically different from that employed in the inversion algorithm. This check provides independent verification of the scheme since the two forward modelling algorithms are prone to different types of numerical error.     

Combination of TOPEX/POSEIDON data with a hydrographic inversion for determination of the oceanic general circulation and its relation to geoid accuracy

A global estimate of the absolute oceanic general circulation from a geostrophic inversion of in situ hydrographic data is tested against and then combined with an estimate obtained from TOPEX/POSEIDON altimetric data and a geoid model computed using the JGM-3 gravity-field solution. Within the quantitative uncertainties of both the hydrographic inversion and the geoid estimate, the two estimates derived by very different methods are consistent. When the in situ inversion is combined with the altimetry/geoid scheme using a recursive inverse procedure, a new solution, fully consistent with both hydrography and altimetry, is found. There is, however, little reduction in the uncertainties of the calculated ocean circulation and its mass and heat fluxes because the best available geoid estimate remains noisy relative to the purely oceano-graphic inferences. The conclusion drawn from this is that the comparatively large errors present in the existing geoid models now limit the ability of satellite altimeter data to improve directly the general ocean circulation models derived from in situ measurements. Because improvements in the geoid could be realized through a dedicated spaceborne gravity recovery mission, the impact of hypothetical much better, future geoid estimates on the circulation uncertainty is also quantified, showing significant hypothetical reductions in the uncertainties of oceanic transport calculations, Full ocean general circulation models could better exploit both existing oceanographic data and future gravity-mission data, but their present use is severely limited by the inability to quantify their error budgets.     

Simultaneous velocity and interface tomography of normal-incidence and wide-aperture seismic traveltime data

A tomographic inversion technique that inverts traveltimes to obtain a model of the subsurface in terms of velocities and interfaces is presented. It uses a combination of refraction, wide-angle reflection and normal-incidence data, it simultaneously inverts for velocities and interface depths, and it is able to quantify the errors and trade-offs in the final model. The technique uses an iterative linearized approach to the non-linear traveltime inversion problem. The subsurface is represented as a set of layers separated by interfaces, across which the velocity may be discontinuous. Within each layer the velocity varies in two dimensions and has a continuous first derivative. Rays are traced in this medium using a technique based on ray perturbation theory, and two-point ray tracing is avoided by interpolating the traveltimes to the receivers from a roughly equidistant fan of rays. The calculated traveltimes are inverted by simultaneously minimizing the misfit between the data and calculated traveltimes, and the roughness of the model. This 'smoothing regularization' stabilizes the solution of the inverse problem. In practice, the first iterations are performed with a high level of smoothing. As the inversion proceeds, the level of smoothing is gradually reduced until the traveltime residual is at the estimated level of noise in the data. At this point, a minimum-feature solution is obtained, which should contain only those features discernible over the noise.
The technique is tested on a synthetic data set, demonstrating its accuracy and stability and also illustrating the desirability of including a large number of different ray types in an inversion.     

Refining seismic structure models by the systematic phase identification of late-arriving groups

We develop a systematic approach to the phase identification of late-arriving groups in 2-D seismic data. Waveforms in the same traveltime branch are grouped, and synthetic traveltimes for all phases are calculated using an initial approximation to the 2-D structure. For each group, we identify the two synthetic phases providing the smallest RMS residuals. If their ratio is less than some predetermined threshold, then the group's phase is ambiguous and both assignments must be tested by traveltime inversion. If there are n unidentified groups, we construct 2 ⁿ phase tables and perform a traveltime inversion on every plausible phase assignment. The phase table that provides the highest value of the posterior probability density is taken as correct, and a 2-D velocity model is constructed from the data. This approach is shown to be effective and efficient on both simulated and real data. In addition, the residuals associated with late-arriving groups provide a means of identifying deficiencies in the initial model.     

Trajectory reconstruction from mobile positioning data using cell-to-cell travel time information

This paper proposes a technique for improving the accuracy of mobile device movement trajectory reconstruction using passive mobile positioning data. The major sources of uncertainty in trajectory reconstruction are imprecise cell shape data and ‘ping-pong’ effects caused by cell handovers. We used a novel technique for improved ‘ping-pong’ effect suppression by compensating for some cell shape distortions based on temporal cell-to-cell transit statistics. The results were evaluated by estimating traffic flow using trajectory reconstruction. The proposed technique improved the accuracy of results compared to ‘ping-pong’ suppression algorithms found in the literature.     

Investigation of a Tertiary maar structure using three-dimensional resistivity imaging

This paper presents the application of the electrical resistivity tomography (ERT) method to the investigation of the Tertiary maar structure of Baruth (Germany) known from previous gravimetric surveys. ERT was applied to support the optimum location for a palaeoclimatological drill hole.
  Special modifications of data acquisition, signal processing and inversion are introduced to adapt the method of ERT to the special requirements for the 3-D investigation of structures with horizontal extensions of 1  km or more. More than 5000 dipole–dipole combinations were recorded at three concentric circular electrode arrangements using stand-alone transient data acquisition systems (RefTek).
  We present a fast approximate imaging technique based on the simultaneous iterative reconstruction technique (SIRT). As the complete calculation of the inverse Frechét matrix is avoided, the algorithm is especially suitable for large data and model spaces, where complete inversion is beyond the limits of available computing hardware. The single-step method is applicable to arbitrary irregular electrode layouts. Synthetic tests show that the imaging procedure reconstructs the main features of the subsurface.
  A low-resistivity body could be interpreted as limnic sediments filling the interior of the Tertiary maar crater. Considering the horizontal resistivity gradient, estimates for the lateral and depth extents of the structure were made. An optimum position for a palaeoclimatological borehole was found, and was in good agreement with the gravimetric minimum.     

Inversion of controlled-source seismic tomography and gravity data with the self-adaptive wavelet parametrization of velocities and interfaces

A self-adaptive automated parametrization approach is suggested for the sequential inversion of controlled-source seismic tomography and gravity data. The velocities and interfaces are parametrized by their Haar wavelet expansion coefficients. Only those coefficients that are well constrained by the data, as measured by the number of rays that cross the corresponding wavelet function support area and their angular coverage, are inverted for, others are set to zero. This approach results in a reasonable distribution of resolution throughout the model even in cases of irregular ray coverage and does overcome the trade-off between different types of model parameters. A modified sequential inversion approach is suggested to join the traveltimes and gravity anomalies inversion. An algorithm is developed that inverts for smooth velocity and density variations inside the seismic layer, the position of its bottom interface as well as for optimal values of the velocity-to-density regression coefficients. The algorithm makes use of direct (diving), reflected and head (critically refracted) wave traveltimes. The algorithm workflow is demonstrated on a synthetic data example.     

Inversion for interface structure using teleseismic traveltime residuals

An inversion method is presented for the reconstruction of interface geometry between two or more crustal layers from teleseismic traveltime residuals. The method is applied to 2-D models consisting of continuous interfaces separating constant-velocity layers. The forward problem of determining ray paths and traveltimes between incident wave fronts below the structure and receivers located on the Earth's surface is solved by an efficient and robust shooting method. A conjugate gradient method is employed to solve the inverse problem of minimizing a least-squares type objective function based on the difference between observed and calculated traveltimes. Teleseismic data do not accurately constrain average vertical structure, so a priori information in the form of layer velocities and average layer thicknesses is required. Synthetic tests show that the method can be used to reconstruct interface geometry accurately, even in the presence of data noise. Tests also show that, if layer velocities and initial interface positions are poorly chosen, lateral structure is still recoverable. The inversion method was applied to previously published teleseismic data recorded by an in-line array of portable seismographs that traversed the northern margin of the Musgrave Block, central Australia. The solution based on interface parametrization is consistent with models given by other studies that used the same data but different methods, most notably the standard tomographic approach that inverts for velocity rather than interface structure.