Estimation Of Cole-Cole Parameters From Time-domain .
Estimation of Cole-Cole parameters from time-domain electromagnetic dataLaurens Beran and Douglas Oldenburg, University of British Columbia.SUMMARYWe present algorithms for the inversion of time-domain electromagnetic data for the parameters of a Cole-Cole model(conductivity, time constant, and chargeability). We applyglobal optimization and nonlinear uncertainty appraisal to aparametric two-layer model, and also develop an underdetermined inversion algorithm. The algorithms are applied to realtime-domain data from Svalbard, Norway.INTRODUCTIONThe occurrence of negative transients in central loop time domain electromagnetic soundings is a well-documented phenomenon in the geophysical literature. While early papersspeculated that this effect could be caused by particular conductivity distributions or magnetic effects, Weidelt (1982) showedthat for a coincident loop system with a step-off primary field,the measured secondary field must be non-negative over nonpolarizable ground, regardless of the subsurface distribution ofconductivity. Observed negative transients can therefore onlybe attributed to polarization effects; that is, the conductivity ofthe ground is time-varying.Subsequent work by Smith (1988) identified “favorable coupling conditions” required for negative transients to be observed in time-domain electromagnetic data:1. The transmitter must be sufficiently close to a polarizable conductor to induce large early-time eddy currents.2. The conductor must be sufficiently polarizable and thereceiver sufficiently close to the conductor that polarization currents can be measured.3. The induced currents must decay sufficiently quicklythat the magnitude of the polarization current can dominate, thereby producing a sign reversal in the measured secondary field.A number of case studies have identified particular geometries which satisfy these favorable coupling conditions. Forexample, Smith and West (1989) consider polarization in a thinconductive layer over a resistive halfspace. In this geometry,strong early-time eddy currents diffuse into the resistive halfspace, where they quickly attenuate. This allows polarizationcurrents in the conductor to dominate the secondary field andproduce a measurable negative transient.To date, forward modelling of observed data remains the extentof data analysis. This is in part because Smith’s coupling conditions are not usually satisfied with commercial time domainsystems, and so there has been little motivation to develop aninversion algorithm for negative transient time-domain data.For example, in systems where receiver and transmitter are notconcentric, polarization effects are more difficult to detect because reversals of the secondary field can occur without polarization. However, negative transients are not uncommonfor helicopter-borne EM systems with central loop configurations (e.g. AeroTEM, VTEM), and this has motivated usto investigate inversion methodologies for processing negativetransients. Routh (1999) developed an inversion algorithm torecover a two-dimensional distribution of polarization parameters from frequency domain resistivity data. In this case, theproblem can be somewhat simplified by estimating a conductivity at each frequency independently. The causal nature oftime-domain data precludes separate inversions of each channel, and so we adopt a functional representation of the conductivity.In this work we investigate methods for inverting negative transients for the parameters of a time-dependent conductivity modelin a one-dimensional layered earth. We begin by presentingthe Cole-Cole model of a time-dependent conductivity, andthen compare overdetermined and underdetermined inversionmethodologies for estimating the parameters of this model.FORWARD MODELLING NEGATIVE TRANSIENTSThe Cole-Cole model is often used to represent a complex conductivity when modelling electromagnetic data. Its frequencydomain form is σ (ω) σ 1 η1 (iωτ)c ,(1)with σ the high-frequency (early time) limit of the conductivity, τ a time constant, and η the chargeability. Both the chargeability and time constant have been found experimentally to bediagnostic of rock type, while time-domain data is quite insensitive to the coefficient c ( Smith (1988)). Following this work,we fix c 1 in our forward modellings and inversions.To compute the time-domain response of a one-dimensionallayered earth with complex conductivity, we use the algorithmdescribed in Farquharson and Oldenburg (1993). The forwardmodelling is computed in the frequency domain at a number offrequencies, with the conductivity at each frequency given byequation 1. The frequency domain solution is then convertedto the time domain using digital filters.Before proceeding to inversion of observed data with this forward model, we investigate the dependence of predicted dataupon the parameters of a polarizable overburden model. Figures 1 and 2 show the dependence of the time of reversal ofthe vertical secondary magnetic field (zero crossing) and theminimum amplitude of the negative portion of the predicteddecay as a function of polarizable overburden thickness andconductivity, respectively. The polarization parameters of theoverburden are fixed at η 0.5, τ 8.5 10 3 and c 1. In
figure 1, the conductivity of the overburden is fixed at 10 1S/m, and in figure 2, the overburden thickness is fixed at 20m.The conductivity of the basement is fixed at σ 10 2 S/m.The source is a step-off waveform transmitted by a wire loopon the ground with radius 100 m.Figure 1: Time of zero crossing and minimum amplitude ofnegative transients as a function of polarizable overburdenthicknessFigure 2: Time of zero crossing and minimum amplitude ofnegative transients as a function of polarizable overburdenconductivity.These figures illustrate the favorable coupling conditions enumerated above. In figure 1, a thicker overburden pushes thetime of polarity reversal later in time. This is because eddycurrents take longer to propagate through a thicker overburden; polarization currents can only dominate the response oncesmoke rings are attenuated in the resistive basement. For thissurvey configuration, the amplitude the minimum of the magnetic field (i.e. the maximum amplitude of the negative transient) in figure 1 shows a distinct minimum. If the polarizable overburden is thin, then there is less polarizable materialpresent and the amplitude of the negative is reduced. As theoverburden thickens, the polarization response grows to a maximum amplitude (minimum value) and is thereafter reduced aseddy currents persist over longer delay times. Variations in theconductivity of the overburden show a similar resonance effect (figure 2). Increasing conductivity initially produces earlier transients as there is more conduction current to chargepolarizable material. However, if the conductivity is too large,conduction currents persist in the overburden and the negativetransient moves later in time. The amplitude of the transientalso shows a nonlinear dependence upon the conductivity ofthe overburden. These examples show how the observation ofa negative transient depends on an optimal coupling of the polarizable body with the sensor. However, as our inversion results in the next sections demonstrate, the nonlinearity of ColeCole model produces strong nonuniquenesses so that multiplemodels can explain the data even when considering a simpletwo-layered earth.OVERDETERMINED INVERSIONIn formulating an inversion procedure for estimating Cole-Coleparameters from negative transient data, we first consider aparametric, overdetermined approach. We choose this approachbased upon Smith’s observation that a thin polarizable layerover a resistive basement is the most common conductivitystructure to produce negative transients. We apply this algorithm to data acquired on Bakaninbreen glacier in Svalbard,Norway. The observed data are shown in figure 4, along withinversion results. Late time data show a characteristic reversal of sign, corresponding to a reversal in the direction of thevertical secondary magnetic field. There is also a second polarity reversal to a positive secondary field at approximately20 ms, though these data are near the estimated noise floor.A one-dimensional layered earth model is consistent with ourFigure 3: Layered structure of Bakaninbreen glacier, Svalbard(from Burgi (2005)).a priori information regarding the structure of the glacier onwhich these data were acquired, and a two-layered model hasbeen used to infer the boundary between clean and sedimentrich ice ( Burgi (2005)). Our model m is parameterized by anoverburden (subscript 1) and basement (subscript 2)m [σ1 , σ2 , τ1 , τ2 , η1 , η2 , t1 ]T(2)with t1 the thickness of the overburden. The nonlinear dependence of the predicted data upon the model parameters discussed in the previous section suggests that the data misfitfunctionφd kWd (dobs F(m))k2(3)may have local minima. In the above expression, the dataweighting matrix Wd quantifies the uncertainty in the observeddata dobs , and F(m) is the data predicted by forward modelling. We found that linearized inversions on synthetic examples using this parameterization converged to local minimaof φd . This problem can be addressed by repeated linearizedinversions from judiciously chosen starting models. Here we
choose instead to minimize the data misfit function with simulated annealing. This is a global optimization algorithm whichuses random perturbations of the model to search model space.Perturbations which decrease the misfit function are always accepted, while perturbations which increase the misfit are accepted probabilistically according to the Metropolis criterion,which is governed by a temperature parameter. This criterionallows the algorithm to escape potentially suboptimal minima.As the temperature parameter is decreased, the model converges to a final model as perturbations are increasingly rejected. Figure 4 shows the convergence of the misfit for accepted models during the simulated annealling search.A straightforward extension of the simulated annealling algorithm, the Metropolis-Hastings sampler, allows us to quantifythe uncertainty in the model parameters. Formally, this algorithm performs numerical integration of the posterior probability density function in N dimensional model space to obtainthe posterior marginal density of the model parameter mi . Ourimplementation of this algorithm follows work by Dosso employing parallel independent samplers to ensure that the posterior density is adequately sampled ( Dosso (2003)). We alsouse models from the simulated annealing inversion to generate perturbations during the sampling process. Figure 4 showssamples from the marginal posterior density functions for theparameters of the parametric model. We see that there are twodistinct peaks of the posterior density function. Surprisingly,both models indicate that the overburden is less conductivethan the basement. In model A, polarizable material is placedat depth, while model B has strong chargeability near the surface. Both models indicate that the overburden thickness isbetween 50 and 60 m, consistent with the cold ice layer shownin figure 3. The ambiguity regarding the location of chargeablematerial is consistent with a previous forward modelling studyof the same data by Burgi (2005).UNDERDETERMINED INVERSIONThe inverse problem can also be tackled with an underdetermined formulation. In this approach we discretize the earthinto a large number of layers and estimate the polarizationparameters (conductivity, time constant, and chargeability) ineach layer. We then minimize the misfit functionφ φd β φm(4)with β a regularization parameter governing the trade-off between data misfit (φd ) and model norm (φm ) components. Asin Farquharson et al. (2003), the model norm constrains themodel via discretized operators W·s (smallness) and W·z (smoothness):φm αsσ kWσs (mσ mσ ,ref )k2 αzσ kWσz (mσ mσ ,ref )k2 αs kWs (mη mη,ref )k2 αz kWz (mη mη,ref )k2ηηηηαsτ kWτs (mτ mτ,ref )k2 αzτ kWτz (mτ mτ,ref )k2 .(5)A priori information regarding model parameters is built intothe inversion via the reference model mref . The coefficientsαs· and αz· balance the contributions of smallness and smoothness components of model parameters to the full model norm.Minimizing the linearized objective function with respect tothe model parameters leads to an overdetermined system ofequations for the model perturbation at each iteration of thealgorithm. This minimization requires the Jacobian matrix ofsensitivities, which we compute analytically in the same manner as Farquharson and Oldenburg (1993).Figure 5 shows an inversion result obtained with this algorithmapplied to the same data considered previously. The fit to thedata is comparable to that obtained by overdetermined inversion. The regularization parameter was “cooled” at each iteration and in figure 5 we have chosen a model from the resulting Tikhonov curve (φd vs. φm ) which adequately balancesthe trade-off between misfit and model norm. The coefficientsα have been selected so that conductivity and chargeabilityhave equal weightings, while the time constant was assigneda weightings three orders of magnitude smaller. Automatedmethods of determining these coefficients have been developedfor similar problems, and this is an avenue for future investigation. In this example, smoothness and smallness terms foreach parameter are weighted equally. This choice of weightings is perhaps responsible for the convergence of the conductivity model to its reference value. In the absence of conductivity structure, the inversion introduces a region of non-zerochargeability down to a depth of 50 m in order to reproduce thenegative portion of the decay. Note also that where the chargeability is zero, the time constant has no effect on the predicteddata, and so we do not display τ for layers where there is negligible chargeability.CONCLUSIONSIn this work we have developed inversion algorithms for recovering one-dimensional models of Cole-Cole parameters fromnegative transient electromagnetic data. Nonlinear overdetermined inversion and uncertainty appraisal of a negative transient sounding acquired on a glacier yielded two models whichpredicted the observed data. One model placed polarizablematerial near the surface, while the second had polarizablematerial at depth. An underdetermined inversion of the samedata placed polarizable material near the surface. All modelsyielded an overburden thickness of approximately 50 m, consistent with prior information regarding the thickness of thecold ice layer in figure 3. Given the ambiguity of the overdetermined inversion result, and the inherent non-uniqueness ofthe underdetermined inversion, we cannot conclude whetherthe source of polarization currents is in this overburden layeror at depth. Soundings with receivers placed outside the loopand horizontal components of the secondary field may helpresolve this ambiguity, and we will investigate this in futurework.ACKNOWLEDGMENTSWe thank Dr. Tavi Murray for providing the data used in thisstudy.
Figure 4: Top row: (l) Misfit as a function of accepted models for simulated annealing inversion of negative transient data. (r)Observed data and data predicted by maximum a posteriori models from uncertainty appraisal. Polarity reversals appear as anincrease in the rate of decay in this plot, with observed negatives beginning at approximately 3 ms and a second reversal at 20 ms.Bottom row: accepted models from nonlinear appraisal. Circle and cross markers indicate the locations of maxima of the posterioridensity function corresponding to models A and B, respectively.Figure 5: Underdetermined inversion of negative transient data. Top left: Tikhonov curve, circle indicates model displayed in thebottom row. Top right: Observed and predicted data. Bottom row: recovered parameters. The reference models for the parameterswere σ 10 2 S/m, τ 10 4 s, and η 0.
APPENDIX ATHE SOURCE OF THE BIBLIOGRAPHY@ARTICLE{Weidelt1982a,author {P. Weidelt},title {Response characteristics of coincident loop transient electromagneticsystems},journal {Geophysics},year {1982},volume {47},pages {1325-1330},}@ARTICLE{Smith1989,author {R. S. Smith and G. F. West},title {Field examples of negative coincident -loop transient electromagneticresponses modeled with polarizable half planes},journal {Geophysics},year {1989},volume {54},pages {1491-1498},}@PHDTHESIS{Smith1988b,author {R. S. Smith},title {A plausible mechanism for generating negative coincident -loop transientelectromagnetic responses},school {University of Toronto},year {1988},}@PHDTHESIS{Routh1999,author {P. S. Routh},title {Electromagnetic coupling in frequency domain induced polarization data.},school {University of British Columbia},year {1999},}@ARTICLE{Farquharson1993,author {C. G. Farquharson and D. W. Oldenburg},title {Inversion of time-domain electromagnetic data for a horizontallylayered Earth},journal {Geophysical Journal International},year {1993},volume {114},pages {433-442},}@ARTICLE{Dosso2003,author {S. Dosso},title {Environmental uncertainty in ocean acoustic source localization},journal {Inverse Problems},year {2003},volume {19},pages {419-431},}
@TECHREPORT{Burgi2005,author {H Burgi},title {Modeling of Induced Polarization Effects in Time-Domain ElectromagneticData from a Glacier in Svalbard, Norway},institution {University of British Columbia},year {2005},}@ARTICLE{Farquharson2003,author {Colin G. Farquharson and Douglas W. Oldenburg and Partha S. Routh},title {Simultaneous 1D inversion of looploop electromagnetic data for magneticsusceptibility and electrical conductivity},journal {Geophysics},year {2003},volume {68},pages {1857-1869},}
REFERENCESBurgi, H., 2005, Modeling of induced polarization effectsin time-domain electromagnetic data from a glacier insvalbard, norway: Technical report, University of BritishColumbia.Dosso, S., 2003, Environmental uncertainty in ocean acousticsource localization: Inverse Problems, 19, 419–431.Farquharson, C. G. and D. W. Oldenburg, 1993, Inversion oftime-domain electromagnetic data for a horizontally layered earth: Geophysical Journal International, 114, 433–442.Farquharson, C. G., D. W. Oldenburg, and P. S. Routh, 2003,Simultaneous 1d inversion of looploop electromagneticdata for magnetic susceptibility and electrical conductivity:Geophysics, 68, 1857–1869.Routh, P. S., 1999, Electromagnetic coupling in frequency domain induced polarization data.: PhD thesis, University ofBritish Columbia.Smith, R. S., 1988, A plausible mechanism for generating negative coincident -loop transient electromagnetic responses:PhD thesis, University of Toronto.Smith, R. S. and G. F. West, 1989, Field examples of negativecoincident -loop transient electromagnetic responses modeled with polarizable half planes: Geophysics, 54, 1491–1498.Weidelt, P., 1982, Response characteristics of coincident looptransient electromagnetic systems: Geophysics, 47, 1325–1330.
the Cole-Cole model of a time-dependent conductivity, and then compare overdetermined and underdetermined inversion . equation 1. The frequency domain solution is then converted to the time domain using digital filters. Before proceeding to inversion of observed data with this for-
t. The log-normal distribution is described by the Cole-Cole a, and the mode of the distribution is the time constant of relaxation [Cole and Cole, 1941]. If the Cole-Cole distribution parameter, a, is unity, then there is a single time constant of relaxation and the Cole
The Cole-Cole (II is a number that is often used to describe the divergence of a measured dielectric dispersion from the ideal dispersion exhibited by a Debye type of dielectric relaxation, and is widely . [27] equation, introduced by the Cole brothers [28] in which an additional parameter, the Cole-Cole (Y, is used to characterise the fact .
the Cole–Cole and PLS models, the latter technique giving more satisfactory results. Keywords On-line biomass monitoring In-situ spectroscopy Scanning capacitance (dielectric) spectroscopy Cole–Cole equation PLS Calibration model robustness Introduction Over the last few decades, the field of biotechnology has
Equation 20 is the prove of equation 1 which relate water saturation to cole cole time; maximum cole cole time and fractal dimension. The capillary pressure can be scaled as logSw (Df 3) logPc constant 21 Where Sw the water saturation, Pc the capillary pressure and
and Cole, 1941) (Cole and Cole, 1941) and GEMTIP model (Zhdanov, 2008). These two models have been used in a number of publications for the interpretation of IP data. The parameters of the conductivity relaxation model can be used for discrimination of the different types of rock formations, which is an important goal in mineral and petroleum
A spreadsheet template for Three Point Estimation is available together with a Worked Example illustrating how the template is used in practice. Estimation Technique 2 - Base and Contingency Estimation Base and Contingency is an alternative estimation technique to Three Point Estimation. It is less
Introduction The EKF has been applied extensively to the field of non-linear estimation. General applicationareasmaybe divided into state-estimation and machine learning. We further di-vide machine learning into parameter estimation and dual estimation. The framework for these areas are briefly re-viewed next. State-estimation
Subclause 1.1 to 1.3 excerpted from ANSI A300 (Part 1) – Pruning 1 ANSI A300 standards 1.1 Scope ANSI A300 standards present performance stan-dards for the care and management of trees, shrubs, and other woody plants. 1.2 Purpose ANSI A300 performance standards are intended for use by federal, state, municipal and private entities
