CN101706587A - Method for extracting induced polarization model parameters prospected by electrical method - Google Patents

Abstract

The invention relates to a method for extracting induced polarization model parameters prospected by electrical methods. The method adopts relative phase spectrum and amplitude spectrum to extract four parameters of a Cole-Cole model. With the effective restriction of the relative phase spectrum to electromagnetic coupling effect, only four parameters of a single Cole-Cole model are needed to be extracted; the extraction method is a combined extraction method which combines the MCMC method in random statistical algorithm with the least square method in deterministic algorithm, and the method has high extraction speed and high extraction accuracy, thus effectively enhancing the work efficiency of induced polarization model parameters, reducing extraction time and effectively restricting interference of electromagnetic effects.

Description

A kind of extracting method of induced polarization model parameters prospected by electrical

Technical field

The present invention relates to a kind of geophysical model parameter extracting method, especially relate to Cole-Cole model parameter extraction method in a kind of resistivity prospecting frequency field induced polarization method.

Background technology

In resistivity prospecting, in order to proofread and correct electromagnetic coupling effect better and to estimate induced polarization anomaly, some scholars have proposed spectral induced polarization method (SIP).This method is to discern polarization body character with complex resistivity, and the evaluation to polarization body is developed into quantitatively from qualitative, has improved greatly and has estimated unusual accuracy.It is in quite wide frequency range (2 ^-8～2 ¹⁰Hz) amplitude and the phase place of measurement complex resistivity in the scope, it is few to have overcome general induced polarization method quantity of information, and can not satisfy the shortcoming that the induced polarization anomaly requirement was proofreaied and correct and estimated to electromagnetic coupling effect, can separate, extract sharp electrical effect and galvanomagnetic effect preferably.Spectral induced polarization method is as the method that directly detects the structure oil-gas possibility, and its prospecting prime cost is low, and the cycle is short, instant effect, and the application in petroleum prospecting and petroleum-gas prediction is extensive day by day.The ultimate principle of this method, field data acquisition instrument and operating technique are increasingly mature.But extract the spectrum parameter exactly and then extract respectively and swash electrical effect and galvanomagnetic effect is still an important content current to be studied.Because the spectrum parameter of spectral induced polarization method acquisition can be the information that two hang-ups (discern and and separate sharp electricity and galvanomagnetic effect and estimate induced polarization anomaly) that solve induced polarization method provide preciousness, thereby this method has caused domestic and international colleague's extensive attention.

The W.H.Pelton of phoenix company etc. thinks that except that swashing electrical effect, the also available Cole-Cole model of the low frequency part of galvanomagnetic effect frequency spectrum is described.So when having sharp electricity and galvanomagnetic effect simultaneously, the complex resistivity frequency characteristic of actual measurement can be expressed as two, three so that a plurality of Cole-Cole model frequency dispersion sums.

${\mathrm{\ρ}}_s\left(\mathrm{i\ω}\right)={\mathrm{\ρ}}_{s0}\{1-m_1[1-\frac{1}{1+{\left({\mathrm{i\ω\τ}}_1\right)}^{c_1}}]-m_2[1-\frac{1}{1+{\left({\mathrm{i\ω\τ}}_2\right)}^{c_2}}]+m_3[1-\frac{1}{1+{\left({\mathrm{i\ω\τ}}_3\right)}^{c_3}}\left]\right\}$

According to the relative size of parameter c and τ, can successfully distinguish the frequency spectrum that swashs electricity and galvanomagnetic effect usually.Yet this method also has its certain defective, in general our fundamental purpose is each parameter that swashs electrical effect in order to try to achieve, but when having electromagnetic coupling effect, we can't ignore its existence, though the model by above-mentioned Cole-Cole frequency dispersion sum can address this problem to a certain extent, but owing to will extract the parameter amount of a times even two to three times more, this undoubtedly can strengthen to a great extent and extract difficulty and reduce the degree of accuracy of extracting, simultaneously, its higher degree of accuracy of data demand of also instrument being surveyed.In addition, this model can only guarantee that the low frequency part of galvanomagnetic effect frequency spectrum can use the Cole-Cole model representation, and is then powerless for the medium-high frequency part.

Although the SIP method can provide the relevant identification of mass efficient to swash information electric and galvanomagnetic effect and evaluation induced polarization anomaly with separating, but try to achieve sharp electrical effect Cole-Cole model parameter with respect to original purpose, it is big that present extraction and evaluation method exist workload, the extraction precision is not enough, accuracy of instrument is had relatively high expectations, can't solve the evaluation problem of galvanomagnetic effect fully.

Summary of the invention

The object of the invention is to provide a kind of extracting parameter amount that reduces, and improves and extracts degree of accuracy, and the extraction of compacting electromagnetic coupling effect swashs the method for electrical model parameters

The purpose that the objective of the invention is to be achieved through the following technical solutions:

One of present inventor once proposed the notion of relative phase spectrum in theory, swashed in the electric relative phase spectrometry in multifrequency, and the relative phase spectrum is defined as:

In the formula:

Expression swashs electric phase spectrum; K represents frequency ratio, is constant, k＞1.If adopt 2n series pseudorandom multi-frequency signal as excitation field source, it is 2,4,8,16 etc. multiple that the value of frequency ratio k can be.This depends on the number and the actual needs of predominant frequency in the pseudorandom multi-frequency signal.

According to the character of Fourier transform and the relation of sharp electrical effect and inductive coupling, prove that the relative phase spectrum has following advantage: can eliminate the measuring error that causes because of the asynchronism(-nization) step between transmitter and the receiver, improve accuracy of observation; The first-order linear correction has been carried out in induction to electromagnetic coupled, can suppress inductive coupling; The observation of relative phase spectrum is measured than frequency dispersion rate higher frequency of operation, and consistent on form with phase spectrum, to equal reaction capacity is arranged unusually.Sharp electric relative phase spectrum and phase spectrum by the field actual measurement show that the relative phase spectrum has good pressing result (referring to Fig. 1) to inductive coupling.

The present invention is based on the Cole-Cole model parameter extraction method of relative phase spectrum and spectral amplitude, it is similar that its extraction principle and spectral amplitude-phase spectrum extracts principle, but rely on its effective compacting to galvanomagnetic effect, directly accurately obtain needed sharp electrical effect model parameter by the extraction to single Cole-Cole model relative phase spectrum and spectral amplitude, single complex resistivity Cole-Cole expression formula is as follows:

$\mathrm{\ρ}\left(\mathrm{i\ω}\right)={\mathrm{\ρ}}_0\{1-m[1-\frac{1}{1+{\left(\mathrm{i\ω\τ}\right)}^c}\left]\right\}---\left(2\right)$

In the formula, ρ (i ω) is a polarization rock complex resistivity when supplying alternating current; ρ ₀Be dc resistivity; M is a polarizability; τ is a time constant; C is a frequency correlation coefficient;

Its spectral amplitude A (ω) is:

$A\left(\mathrm{\ω}\right)={\{{\left[\mathrm{Re\ρ}\left(\mathrm{i\ω}\right)\right]}^2+{\left[\mathrm{Im\ρ}\left(\mathrm{i\ω}\right)\right]}^2\}}^{1/2}$

$={\mathrm{\&rho;}}_0{[(1-m+\frac{\mathrm{mR}}{R^2+I^2})+\left(\frac{\mathrm{<figure-callout\; id="2"\; label="mI\; R"\; filenames="H2009102270095E0000022.png"\; state="\{\{state\}\}">mI</figure-callout>}}{R^{}}\right)]}^{1/2}---\left(3\right)$

$={\mathrm{\&rho;}}_0{\left[\frac{1+2(1-m){\left(\mathrm{\&omega;\&tau;}\right)}^c\cos \frac{\mathrm{\&pi;<figure-callout\; id="2"\; label="c"\; filenames="H2009102270095E0000022.png"\; state="\{\{state\}\}">c</figure-callout>}}{2}+{(1-m)}^2{\left(\mathrm{\&omega;\&tau;}\right)}^{2\mathrm{<figure-callout\; id="1"\; label="c"\; filenames="H2009102270095E0000011.png"\; state="\{\{state\}\}">c</figure-callout>}}}{1+2{\left(\mathrm{\&omega;\&tau;}\right)}^c\cos \frac{\mathrm{\&pi;<figure-callout\; id="2"\; label="c"\; filenames="H2009102270095E0000022.png"\; state="\{\{state\}\}">c</figure-callout>}}{2}+{\left(\mathrm{\&omega;\&tau;}\right)}^{2c}}\right]}^{1/2}$

Its phase spectrum

For:

According to the computing formula (1) of formula relative phase spectrum definition, try to achieve the relative phase spectrum; By relative phase spectrum definition, it is as follows at each Parameter Extraction formula of Cole-Cole to get the relative phase spectrum:

In described formula (3), (4), (5), (6), (7):

$\left.\begin{array}{c}{\mathrm{<figure-callout\; id="1"\; label="R"\; filenames="H2009102270095E0000011.png"\; state="\{\{state\}\}">R</figure-callout>}}_1=1+{\left(\mathrm{\&omega;\&tau;}\right)}^c\cos \frac{\mathrm{\&pi;c}}{2}\\ {\mathrm{<figure-callout\; id="1"\; label="I"\; filenames="H2009102270095E0000011.png"\; state="\{\{state\}\}">I</figure-callout>}}_1={\left(\mathrm{\&omega;\&tau;}\right)}^c\sin \frac{\mathrm{\&pi;c}}{2}\end{array}\right\}$ $\left.\begin{array}{c}{\mathrm{<figure-callout\; id="2"\; label="R"\; filenames="H2009102270095E0000022.png"\; state="\{\{state\}\}">R</figure-callout>}}_2=1+{\left(\mathrm{k\&omega;\&tau;}\right)}^c\cos \frac{\mathrm{\&pi;c}}{2}\\ {\mathrm{<figure-callout\; id="2"\; label="I"\; filenames="H2009102270095E0000022.png"\; state="\{\{state\}\}">I</figure-callout>}}_2={\left(\mathrm{k\&omega;\&tau;}\right)}^c\sin \frac{\mathrm{\&pi;<figure-callout\; id="2"\; label="c"\; filenames="H2009102270095E0000022.png"\; state="\{\{state\}\}">c</figure-callout>}}{2}\end{array}\right\}$

${\mathrm{<figure-callout\; id="1"\; label="R"\; filenames="H2009102270095E0000011.png"\; state="\{\{state\}\}">R</figure-callout>}}_1^2+{\mathrm{<figure-callout\; id="1"\; label="I"\; filenames="H2009102270095E0000011.png"\; state="\{\{state\}\}">I</figure-callout>}}_1^2=1+2{\left(\mathrm{\&omega;\&tau;}\right)}^c\cos \frac{\mathrm{\&pi;<figure-callout\; id="2"\; label="c"\; filenames="H2009102270095E0000022.png"\; state="\{\{state\}\}">c</figure-callout>}}{2}+{\left(\mathrm{\&omega;\&tau;}\right)}^{2\mathrm{<figure-callout\; id="2"\; label="c\; R"\; filenames="H2009102270095E0000022.png"\; state="\{\{state\}\}">c</figure-callout>}}$ $R_{}^{}$ ${}_2^2+{\mathrm{<figure-callout\; id="2"\; label="I"\; filenames="H2009102270095E0000022.png"\; state="\{\{state\}\}">I</figure-callout>}}_2^2=1+2{\left(\mathrm{k\&omega;\&tau;}\right)}^c\cos \frac{\mathrm{\&pi;<figure-callout\; id="2"\; label="c"\; filenames="H2009102270095E0000022.png"\; state="\{\{state\}\}">c</figure-callout>}}{2}+{\left(\mathrm{k\&omega;\&tau;}\right)}^{2c}$

${\mathrm{<figure-callout\; id="1"\; label="I"\; filenames="H2009102270095E0000011.png"\; state="\{\{state\}\}">I</figure-callout>}}_1^2-{\mathrm{<figure-callout\; id="1"\; label="R"\; filenames="H2009102270095E0000011.png"\; state="\{\{state\}\}">R</figure-callout>}}_1^2={\left(\mathrm{\&omega;\&tau;}\right)}^{2c}({\mathrm{<figure-callout\; id="2"\; label="sin"\; filenames="H2009102270095E0000022.png"\; state="\{\{state\}\}">sin</figure-callout>}}^2\frac{\mathrm{\&pi;c}}{2}-{\mathrm{<figure-callout\; id="2"\; label="cos"\; filenames="H2009102270095E0000022.png"\; state="\{\{state\}\}">cos</figure-callout>}}^2\frac{\mathrm{\&pi;c}}{2})-1-2{\left(\mathrm{\&omega;\&tau;}\right)}^c\cos \frac{\mathrm{\&pi;<figure-callout\; id="2"\; label="c"\; filenames="H2009102270095E0000022.png"\; state="\{\{state\}\}">c</figure-callout>}}{2}$

${\mathrm{<figure-callout\; id="2"\; label="I"\; filenames="H2009102270095E0000022.png"\; state="\{\{state\}\}">I</figure-callout>}}_2^2-{\mathrm{<figure-callout\; id="2"\; label="R"\; filenames="H2009102270095E0000022.png"\; state="\{\{state\}\}">R</figure-callout>}}_2^2={\left(\mathrm{k\&omega;\&tau;}\right)}^{2c}({\sin }^2\frac{\mathrm{\&pi;c}}{2}-{\mathrm{<figure-callout\; id="2"\; label="cos"\; filenames="H2009102270095E0000022.png"\; state="\{\{state\}\}">cos</figure-callout>}}^2\frac{\mathrm{\&pi;c}}{2})-1-2{\left(\mathrm{k\&omega;\&tau;}\right)}^c\cos \frac{\mathrm{\&pi;<figure-callout\; id="2"\; label="c"\; filenames="H2009102270095E0000022.png"\; state="\{\{state\}\}">c</figure-callout>}}{2}$

$\frac{\&PartialD;{\mathrm{<figure-callout\; id="1"\; label="R"\; filenames="H2009102270095E0000011.png"\; state="\{\{state\}\}">R</figure-callout>}}_1}{\&PartialD;c}={\left(\mathrm{\&omega;\&tau;}\right)}^c[\ln \left(\mathrm{\&omega;\&tau;}\right)\&CenterDot;\cos \frac{\mathrm{\&pi;c}}{2}-\frac{\mathrm{\&pi;}}{2}\sin \frac{\mathrm{\&pi;c}}{2}]$

$\frac{\&PartialD;{\mathrm{<figure-callout\; id="2"\; label="R"\; filenames="H2009102270095E0000022.png"\; state="\{\{state\}\}">R</figure-callout>}}_2}{\&PartialD;c}={\left(\mathrm{k\&omega;\&tau;}\right)}^c[\ln \left(\mathrm{k\&omega;\&tau;}\right)\&CenterDot;\cos \frac{\mathrm{\&pi;c}}{2}-\frac{\mathrm{\&pi;}}{2}\sin \frac{\mathrm{\&pi;c}}{2}]$

$\frac{\&PartialD;I_1}{\&PartialD;c}={\left(\mathrm{\&omega;\&tau;}\right)}^c[\ln \left(\mathrm{\&omega;\&tau;}\right)\&CenterDot;\sin \frac{\mathrm{\&pi;<figure-callout\; id="2"\; label="c"\; filenames="H2009102270095E0000022.png"\; state="\{\{state\}\}">c</figure-callout>}}{2}+\frac{\mathrm{\&pi;}}{2}\cos \frac{\mathrm{\&pi;c}}{2}]$

$\frac{\&PartialD;{\mathrm{<figure-callout\; id="2"\; label="I"\; filenames="H2009102270095E0000022.png"\; state="\{\{state\}\}">I</figure-callout>}}_2}{\&PartialD;c}={\left(\mathrm{k\&omega;\&tau;}\right)}^c[\ln \left(\mathrm{k\&omega;\&tau;}\right)\&CenterDot;\sin \frac{\mathrm{\&pi;<figure-callout\; id="2"\; label="c"\; filenames="H2009102270095E0000022.png"\; state="\{\{state\}\}">c</figure-callout>}}{2}+\frac{\mathrm{\&pi;}}{2}\cos \frac{\mathrm{\&pi;c}}{2}]$

$\frac{{\&PartialD;\mathrm{<figure-callout\; id="1"\; label="R"\; filenames="H2009102270095E0000011.png"\; state="\{\{state\}\}">R</figure-callout>}}_1}{\&PartialD;\mathrm{\&tau;}}=\frac{c}{\mathrm{\&tau;}}{\left(\mathrm{\&omega;\&tau;}\right)}^c\cos \frac{\mathrm{\&pi;<figure-callout\; id="2"\; label="c"\; filenames="H2009102270095E0000022.png"\; state="\{\{state\}\}">c</figure-callout>}}{2}$ $\frac{{\&PartialD;\mathrm{<figure-callout\; id="2"\; label="R"\; filenames="H2009102270095E0000022.png"\; state="\{\{state\}\}">R</figure-callout>}}_2}{\&PartialD;\mathrm{\&tau;}}=\frac{c}{\mathrm{\&tau;}}{\left(\mathrm{k\&omega;\&tau;}\right)}^c\cos \frac{\mathrm{\&pi;<figure-callout\; id="2"\; label="c"\; filenames="H2009102270095E0000022.png"\; state="\{\{state\}\}">c</figure-callout>}}{2}$

$\frac{{\&PartialD;I}_1}{\&PartialD;\mathrm{\&tau;}}=\frac{c}{\mathrm{\&tau;}}{\left(\mathrm{\&omega;\&tau;}\right)}^c\sin \frac{\mathrm{\&pi;<figure-callout\; id="2"\; label="c"\; filenames="H2009102270095E0000022.png"\; state="\{\{state\}\}">c</figure-callout>}}{2}$ $\frac{{\&PartialD;I}_2}{\&PartialD;\mathrm{\&tau;}}=\frac{c}{\mathrm{\&tau;}}{\left(\mathrm{k\&omega;\&tau;}\right)}^c\sin \frac{\mathrm{\&pi;<figure-callout\; id="2"\; label="c"\; filenames="H2009102270095E0000022.png"\; state="\{\{state\}\}">c</figure-callout>}}{2}$

Because the relative phase spectrum only includes 3 unknown parameters, therefore a remaining unknown parameter ρ ₀Need to rely on spectral amplitude to extract, because the two is linear, it is as follows to draw its extraction fundamental formular according to formula (3):

$\frac{\&PartialD;A\left(\mathrm{\&omega;}\right)}{{\&PartialD;\mathrm{\&rho;}}_0}={\left[\frac{1+2(1-m){\left(\mathrm{\&omega;\&tau;}\right)}^c\cos \frac{\mathrm{\&pi;<figure-callout\; id="2"\; label="c"\; filenames="H2009102270095E0000022.png"\; state="\{\{state\}\}">c</figure-callout>}}{2}+{(1-m)}^2{\left(\mathrm{\&omega;\&tau;}\right)}^{2\mathrm{<figure-callout\; id="1"\; label="c"\; filenames="H2009102270095E0000011.png"\; state="\{\{state\}\}">c</figure-callout>}}}{1+2{\left(\mathrm{\&omega;\&tau;}\right)}^c\cos \frac{\mathrm{\&pi;<figure-callout\; id="2"\; label="c"\; filenames="H2009102270095E0000022.png"\; state="\{\{state\}\}">c</figure-callout>}}{2}+{\left(\mathrm{\&omega;\&tau;}\right)}^{2c}}\right]}^{1/2}---\left(9\right)$

On the other hand, in order to improve the extraction effect of relative phase spectrum and spectral amplitude, we need seek the extracting method that is appropriate to the relative phase spectrum.The technical solution adopted in the present invention is to unite with the least square method in the deterministic algorithm in conjunction with random statistical algorithm MCMC method (the special calot's method of Marko Fumeng) to extract relative phase spectrum and spectral amplitude.The MCMC method that is wherein adopted is used the Metropolis-Hasting algorithm, produces and analyzes the Metropolis-Hasting Markov chain and finally extracted the result by polynary Gaussian distribution method.

Whole technical proposal may further comprise the steps:

1) data that go out relative phase spectrum and spectral amplitude theoretical value with giving Cole-Cole Model Calculation or directly use the measured value that recorded in the instrument to extract as model; Selected Cole-Cole model parameter to be asked is set up the theory relation of waiting to ask between model parameter and relative phase spectrum and spectral amplitude observation data according to extraction model formula (5), (6), (7), (9); According to the actual requirement and the prior imformation of the problem of extraction, select the limits of error between model theory value and the observed reading, be limited to the error of measurement data on the error;

2) use decisive algorithm (traditional least square method): setting up the least square model on original relative phase spectrum and spectral amplitude model based is least square and model, and set by step 1) given standard provides relevant parameter setting in, adopt relative phase spectrum and spectral amplitude data given in the step 1), appoint according to the scope of waiting to ask parameter in the step 3) and to establish one group of parameter initial value, extract fundamental formular according to the relative phase spectrum and calculate near the locally optimal solution of model initial value with least square method;

The method to set up of described parameter initial value is: the resistivity initial value is the resistivity of spectral amplitude low-limit frequency correspondence in the measurement data, the frequency correlation coefficient initial value is 0.2-0.5, the polarizability initial value is that the time constant initial value is the inverse of the frequency of relative phase polarographic maximum correspondence with the maximum frequency dispersion rate of spectral amplitude definition.

3) use MCMC method (randomness algorithm): the initial value that the model optimum solution of being calculated employing step 2) is extracted as the MCMC method, when being provided with, MCMC method parameter provides parameter value priori scope to be measured according to the actual conditions in image data or the application process, and in computing machine, generate the multiple parameter model of Bayes at random by the foundation of priori scope, produce and analyze the Metropolis-Hasting Markov chain with the MCMC method by the Gaussian distribution method then, thereby try to achieve the overall edge posterior probability density function of each unknown parameter in effective range; According in the step 1) the given standard accepted check generation model and calculate the extraction precision;

Described parameter value priori scope to be measured is: resistivity priori scope is a complex resistivity spectral amplitude variation range, frequency correlation coefficient priori value range is 0.001-0.99, polarizability priori scope is 0.001-0.99, and time constant priori scope is 0.001-10000s.

4) testing model: the gained precision does not reach the regulation requirement if extract as a result, and the overall probability density function that extract with MCMC method in step 3) this moment asks the parameter initial value still to carry out step 3) as waiting, so repeats to meet the demands up to extracting precision; If extracting precision reaches requirement and then uses step 2) in the overall probability density function that extracts with the MCMC method of this moment of least square method be initial value, thereby extract the globally optimal solution that obtains on the practical significance again;

5) result to be extracted stable after, carry out match to data and the result analyzed extracting curve and institute at last.

The present invention extracts 4 parameters of Cole-Cole model by relative phase spectrum and spectral amplitude, rely on of the effective compacting of relative phase spectrum to electromagnetic coupling effect, only need to extract 4 parameters of single Cole-Cole model, and the Instrument observation precision of relative phase spectrum is higher, so better on the more former basis of extraction effect to the extraction of a plurality of Cole-Cole model sums, can obtain parameter value accurately substantially.Extracting on the concrete grammar, we once adopted least square method and genetic algorithm to carry out the extraction of relative phase spectrum separately, but since spectrum induced polarization extract in initial value choose and have considerable influence (its convergence is stable inadequately) and genetic algorithm consuming time more least square method, both guaranteed by MCMC method and least square method associating extracting method that not choosing of initial value can impact and obtain a large amount of extraction uncertainty information to extracting the result, also retained simultaneously least square method extraction rate advantage faster, by our experimental results show that, this associating extracting method generally just can get stabilization result to the end through two to three extractions.

Description of drawings

Fig. 1 is phase spectrum (curve 1) and relative phase spectrum (curve 2) curve of open-air actual measurement.

Fig. 2 is that MCMC method of the present invention and least square method are united extraction workflow synoptic diagram.

Fig. 3 is that spectral amplitude is at ρ ₀=25; τ=100; C=0.25; The m=0.5 theoretical value is extracted matched curve.

Fig. 4 composes at ρ for relative phase ₀=25; τ=100; C=0.25; The m=0.5 theoretical value is extracted matched curve.

Fig. 5 is that spectral amplitude is at ρ ₀=25; τ=100; C=0.25; When adding 1% relative noise, the m=0.5 theoretical value extracts matched curve.

Fig. 6 composes at ρ for relative phase ₀=25; τ=100; C=0.25; When adding 1% relative noise, the m=0.5 theoretical value extracts matched curve.

Fig. 7-9 is the extraction matched curve of relative phase spectrum when considering electromagnetic coupling effect.

Figure 10 extracts matched curve for actual measurement pseudorandom 7 frequency wave datum in the relative phase spectrum of exception.

Embodiment

The present invention will be further described below in conjunction with drawings and Examples.

Introduce at first in detail with the extraction flow process of MCMC method to relative phase spectrum and spectral amplitude, its concrete implementation step is as follows:

1. the theoretical value of calculating them according to 20 given class frequency values with the existing Cole-Cole model relative phase spectrum and the formula of spectral amplitude is organized data with the calculated value behind the relative noise of adding as y, 20 groups the frequency values (f=2 of giving ^k, k=-12 ,-11 ,-10 ..., 6,7) and as x group data; If fruit is used 7 given frequency ripple measured data values, then according to the given frequency (f=2 of instrument ^k, k=-2 ,-1,0,1,2,3) and as x group data;

2. relative phase spectrum and spectral amplitude pattern function are provided, in order to accelerate extraction rate, obtain its least square and function according to model, appoint and to establish an initial value in the parameter value scope, utilize the least square method algorithm that provides previously to obtain the one local minimum and begin extraction initial value as the MCMC method; The extraction process step of least square method (Gauss-Newton method) is as follows:

1) given initial value point b ⁽⁰⁾(the initial value vector that contains given unknown parameter) allows to extract error ε＞0, makes iterations k=0;

2) calculate Cole-Cole model component function value f _i(b ^(k)), i=1,2 ..., M gets vectorial

f ^(k)＝[f ₁(b ^(k))，f ₂(b ^(k))，…，f _M(b ^(k))] ^T，

Calculate first order derivative again: $a_{\mathrm{ij}}={\&PartialD;f}_i\left(b^{\left(k\right)}\right)/{\&PartialD;b}_j(i=\mathrm{1,2},...,M;j=\mathrm{1,2},...,N)$

Get MxN rank matrix: A _k={ a _Ij} _MxN

3) group of solving an equation: Get Gauss-newton's direction p ^(k)

4) from b ^(k)Set out, along direction p ^(k)Do linear search, determine optimal step size t ^(k)

5) if || b ^(k+1)-b ^(k)||＜ε, then stop to calculate, get b ^*=b ^(k+1), otherwise k=k+1 returned for the 2nd step in addition;

Begin relevant parameter setting afterwards, each parameter range is set, earlier as ρ ₀Scope is roughly 0-100, and the scope of τ is set to 0-1000, and the c scope is set to 0-1, and the scope of m is set to 0-1.The extraction precision then is set, extracts calculation times and extract sampling parameters or the like;

3. just can directly begin to extract with the MCMC method after setting up parameter, be example in conjunction with the concrete workflow synoptic diagram of Fig. 2 with the relative phase spectrum, establish 3 unknown parameters and be respectively m, and b, c, it is as follows that the MCMC method is extracted idiographic flow:

(1) tries to achieve the relative phase spectrum by Gaussian distribution Likelihood equation:

μ is its relative error variable

Ask the conditional probability distribution function (is example with m) of its each unknown parameter then:

Ind (m ∈ D _m) expression extracts the span at m place;

(2) initialize is given m, b, and c and μ can write a Chinese character in simplified form into m respectively ⁽⁰⁾, b ⁽⁰⁾, c ⁽⁰⁾And μ ⁽⁰⁾And allow k=1;

(3) sampling and given b from f (m|.) ^(k-1), c ^(k-1)And μ ^(k-1)And it is abbreviated as m ^(k)

(4) sampling and given m from f (b|.) ^(k), c ^(k-1)And μ ^(k-1)And it is abbreviated as b ^(k)

(5) sampling and given m from f (c|.) ^(k), b ^(k)And μ ^(k-1)And it is abbreviated as c ^(k)

(6) from f (μ | .) the sampling and given m ^(k), b ^(k)And c ^(k)And it is abbreviated as μ ^(k)

(7) allow k=k+1, if k＞T, wherein T is the maximum iterations that allows, and then stops; Otherwise, execution in step (2);

(8) the last result that extracts of output;

Program can provide each parameter value after the extraction, because can't determine whether fully that parameter value is near actual value, therefore the parameter value after extracting need be extracted again as the initial value in the 2nd step once again, if have parameter value to exceed its given scope after extracting, then when extracting next time, still use that initial parameter initial value.So repeat up to extract the result stable after or its standard deviation that calculates can stop extraction within error range the time.Because Markov chain is that overall probability is chosen random point, therefore what obtain is the posterior probability density function of the overall situation, chooses regardless of initial value, only can influence that its extraction time produces and can not influence it and extract the result, in general, 2-3 extraction can reach end product;

4. after extracting end, draw its model points and also calculate some parameters accordingly, and analyze corresponding results, extraction effect is estimated with the matched curve of extracting the result;

The present invention adopts the MCMC method to 4 parameter ρ in the Cole-Cole model ₀, τ, c, m extracts calculating, ρ ₀, τ, c, that the theoretical value of m adopts is Luo Yanzhong, the reference value that Zhang Guiqing provides, and frequency ratio k value all is 2 in the relative phase spectrum, because not choosing of parameter initial value can change the extraction net result, only can influence the time of its extraction, so our concrete implementation step is as follows:

1. spectral amplitude and relative phase spectral theory calculated value are extracted, check that can the MCMC method reach our required extraction accuracy requirement; So need obtain another parameter ρ by extracting spectral amplitude because the relative phase spectrum only contains 3 parameters ₀Value; Fig. 3 is that spectral amplitude is at theoretical value ρ ₀=25; τ=100; C=0.25; The matched curve of extracting on the m=0.5 basis sees that by knowing among the figure extraction degree of fitting is very high, and the extraction result that the MCMC method provides is: ρ ₀=24.9997; τ=99.9822; C=0.25001; M=0.49999, the extraction error is 7.7543e-009.Fig. 4 is the matched curve that the relative phase spectrum is extracted on the theoretical value basis, and its degree of fitting is equally very high, and the extraction result that its MCMC method provides is: τ=100.19; C=0.250032; M=0.499985, the extraction error is 1.0425e-008.According to above-mentioned experimental result, the extraction precision of MCMC method reaches necessary requirement;

2. the theoretical value of spectral amplitude and relative phase spectrum is added that 1% relative level of noise extracts, the value of noise is with reference to interrelated data relatively; Carrying out this step experiment is in order to detect relative phase spectrum and spectral amplitude the extraction precision of using the MCMC method under the noise jamming to be arranged, and theoretical value is still continued to use the value in the 1st step; Fig. 5 has reflected the matched curve that the spectral amplitude theoretical value is extracted when the relative noise that adds 1%, curve fitting is better, and its extraction result is: ρ ₀=24.2889; τ=56.6953; C=0.252691; M=0.48851, extracting error is 0.3482.From extracting on curve and the result, except time constant, other parameter extraction all compares near actual value, and this also conforms to the relatively poor fact of spectral amplitude extraction time constant effect; Fig. 6 has reflected the extraction matched curve when the relative noise that adds 1% of relative phase spectral theory value, and the curve fitting degree is still higher, and its extraction result is: τ=96.1673; C=0.251792; M=0.498894, the extraction error is 3.8914e-006, very near actual value, it can well remedy the weak point of spectral amplitude to visible relative phase spectrum to three Parameter Extraction results;

3. according to relative phase spectrum noted earlier electromagnetic coupling effect there is suppression, we have added the value of that part of Cole-Cole model of electromagnetic coupled in theoretical value is calculated, extracting with single Cole-Cole model relative phase spectrum then, by extracting the compacting electromagnetic coupling effect effect that the result confirms that relative phase is composed. the electromagnetic coupled Cole-Cole model parameter value that we choose also is a reference value: τ=1; C=0.98; M=0.01. Fig. 7 is τ=100 for getting theoretical value; C=0.25; Extraction matched curve during m=0.5, its MCMC method are extracted the result and are: τ=71.203; C=0.252177; M=0.510057, extracting error is τ=10 for 1.0889e-004. Fig. 8 for getting theoretical value; C=0.4; Extraction matched curve during m=0.5, its MCMC method are extracted the result and are: τ=8.87723; C=0.404417; M=0.509939, extracting error is τ=10 for 6.5350e-005. Fig. 9 for getting theoretical value; C=0.6; Extraction matched curve during m=0.9, its MCMC method are extracted the result and are: τ=9.5103; C=0.601493; M=0.910012, extracting error is the extraction result of 2.3529e-004. from above several groups of data, the relative phase spectrum is extracted and can be suppressed the influence of galvanomagnetic effect really, improves the extraction precision that swashs electrical effect;

4. extract with pseudorandom multi-frequency induced polarization instrument one group of 7 measured frequency ripple measured data in flume experiment, because frequency band range is smaller from 0.25Hz-8Hz, so can only measure a part of curve of relative phase spectrum.Be placed with copper coin in the tank, the extraction matched curve of the relative phase spectrum of carrying out according to measured data as shown in figure 10, it extracts the result and is: τ=7.29688e+010; C=0.151606; M=0.916147, extracting error is 0.000014, because in the low-resistivity exception, the time constant value has been difficult to accurate Calculation, therefore above-mentioned data can only reflect that low-resistivity is that copper coin is unusual unusually, c and m approach actual value.

Above explanation mode only is used to illustrate the present invention; and be not limiting the scope of the invention; the those of ordinary skill in relevant technologies field; under the situation that does not break away from the spirit and scope of the present invention; can also make various variations and distortion, all technical schemes that are equal to also belong to category of the present invention.

Claims (1)

1. the extracting method of an induced polarization model parameters prospected by electrical, it is characterized in that, utilize relative phase spectrum and spectral amplitude to extract the Cole-Cole model parameter, obtain required sharp electrical effect model parameter by extraction to single Cole-Cole model relative phase spectrum and spectral amplitude;

The method of utilizing relative phase spectrum and spectral amplitude to extract the Cole-Cole model parameter is: in conjunction with the least square method in random statistical algorithm MCMC method and the deterministic algorithm, extract the Cole-Cole model parameter from relative phase spectrum and spectral amplitude; The MCMC method that is wherein adopted is used the Metropolis-Hasting algorithm, produces and analyzes the Metropolis-Hasting Markov chain and finally extracted the result by polynary Gaussian distribution method;

The step that obtains required sharp electrical effect model parameter by the extraction to single Cole-Cole model relative phase spectrum and spectral amplitude is:

1) data that go out relative phase spectrum and spectral amplitude theoretical value with giving Cole-Cole Model Calculation or directly use the measured value that recorded in the instrument to extract as model; Selected Cole-Cole model parameter to be asked, set up the theory relation of waiting to ask between model parameter and relative phase spectrum and spectral amplitude observation data according to following extraction model formula (5), (6), (7), (9):

In described formula (5), (6), (7):

$\left.\begin{array}{c}{R}_1=1+{\left(\mathrm{\&omega;\&tau;}\right)}^c\cos \frac{\mathrm{\&pi;c}}{2}\\ I_1={\left(\mathrm{\&omega;\&tau;}\right)}^c\sin \frac{\mathrm{\&pi;c}}{2}\end{array}\right\}$ $\left.\begin{array}{c}{R}_2=1+{\left(\mathrm{k\&omega;\&tau;}\right)}^c\cos \frac{\mathrm{\&pi;c}}{2}\\ I_2={\left(\mathrm{k\&omega;\&tau;}\right)}^c\sin \frac{\mathrm{\&pi;c}}{2}\end{array}\right\}$

$\begin{array}{cc}{R}_1^2+I_1^2=1+2{\left(\mathrm{\&omega;\&tau;}\right)}^c\cos \frac{\mathrm{\&pi;c}}{2}+{\left(\mathrm{\&omega;\&tau;}\right)}^{2c}& R_2^2+I_2^2=1+2{\left(\mathrm{k\&omega;\&tau;}\right)}^c\cos \frac{\mathrm{\&pi;c}}{2}+{\left(\mathrm{k\&omega;\&tau;}\right)}^{2c}\end{array}$

$I_1^2-R_1^2={\left(\mathrm{\&omega;\&tau;}\right)}^{2c}=({\sin }^2\frac{\mathrm{\&pi;c}}{2}-c{\mathrm{os}}^2\frac{\mathrm{\&pi;c}}{2})-1-2{\left(\mathrm{\&omega;\&tau;}\right)}^c\cos \frac{\mathrm{\&pi;c}}{2}$

$I_2^2-R_2^2={\left(\mathrm{k\&omega;\&tau;}\right)}^{2c}({\sin }^2\frac{\mathrm{\&pi;c}}{2}-{\cos }^2\frac{\mathrm{\&pi;c}}{2})-1-2{\left(\mathrm{k\&omega;\&tau;}\right)}^c\cos \frac{\mathrm{\&pi;c}}{2}$

$\frac{\&PartialD;R_1}{\&PartialD;c}={\left(\mathrm{\&omega;\&tau;}\right)}^c[\ln \left(\mathrm{\&omega;\&tau;}\right)\&CenterDot;\cos \frac{\mathrm{\&pi;c}}{2}-\frac{\mathrm{\&pi;}}{2}\sin \frac{\mathrm{\&pi;c}}{2}]$

$\frac{\&PartialD;R_2}{\&PartialD;c}={\left(\mathrm{k\&omega;\&tau;}\right)}^c[\ln \left(\mathrm{k\&omega;\&tau;}\right)\&CenterDot;\cos \frac{\mathrm{\&pi;c}}{2}-\frac{\mathrm{\&pi;}}{2}\sin \frac{\mathrm{\&pi;c}}{2}]$

$\frac{\&PartialD;I_1}{\&PartialD;c}={\left(\mathrm{\&omega;\&tau;}\right)}^c[\ln \left(\mathrm{\&omega;\&tau;}\right)\&CenterDot;\sin \frac{\mathrm{\&pi;c}}{2}+\frac{\mathrm{\&pi;}}{2}\cos \frac{\mathrm{\&pi;c}}{2}]$

$\frac{\&PartialD;I_2}{\&PartialD;c}={\left(\mathrm{k\&omega;\&tau;}\right)}^c[\ln \left(\mathrm{k\&omega;\&tau;}\right)\&CenterDot;\sin \frac{\mathrm{\&pi;c}}{2}+\frac{\mathrm{\&pi;}}{2}\cos \frac{\mathrm{\&pi;c}}{2}]$

$\frac{\&PartialD;R_1}{\&PartialD;\mathrm{\&tau;}}=\frac{c}{\mathrm{\&tau;}}{\left(\mathrm{\&omega;\&tau;}\right)}^c\cos \frac{\mathrm{\&pi;c}}{2}$ $\frac{\&PartialD;R_2}{\&PartialD;\mathrm{\&tau;}}=\frac{c}{\mathrm{\&tau;}}{\left(\mathrm{k\&omega;\&tau;}\right)}^c\cos \frac{\mathrm{\&pi;c}}{2}$

$\frac{\&PartialD;I_1}{\&PartialD;\mathrm{\&tau;}}=\frac{c}{\mathrm{\&tau;}}{\left(\mathrm{\&omega;\&tau;}\right)}^c\sin \frac{\mathrm{\&pi;c}}{2}$ $\frac{\&PartialD;I_2}{\&PartialD;\mathrm{\&tau;}}=\frac{c}{\mathrm{\&tau;}}{\left(\mathrm{k\&omega;\&tau;}\right)}^c\sin \frac{\mathrm{\&pi;c}}{2}$

$\frac{\&PartialD;A\left(\mathrm{\&omega;}\right)}{\&PartialD;{\mathrm{\&rho;}}_0}={\left[\frac{1+2(1-m){\left(\mathrm{\&omega;\&tau;}\right)}^c\cos \frac{\mathrm{\&pi;c}}{2}+{(1-m)}^2{\left(\mathrm{\&omega;\&tau;}\right)}^{2c}}{1+2{\left(\mathrm{\&omega;\&tau;}\right)}^c\cos \frac{\mathrm{\&pi;c}}{2}+{\left(\mathrm{\&omega;\&tau;}\right)}^{2c}}\right]}^{1/2}---\left(9\right)$

According to the actual requirement and the prior imformation of the problem of extraction, select the limits of error between model theory value and the observed reading, be limited to the error of measurement data on the error;

2) using traditional least square method to set up the least square model on original relative phase spectrum and spectral amplitude model based is least square and model, and set by step 1) given standard provides relevant parameter setting in, adopt relative phase spectrum and spectral amplitude data given in the step 1), appoint according to the scope of waiting to ask parameter in the step 3) and to establish one group of parameter initial value, extract fundamental formular according to the relative phase spectrum and calculate near the locally optimal solution of model initial value with least square method;

The method to set up of described parameter initial value is: the resistivity initial value is the resistivity of spectral amplitude low-limit frequency correspondence in the measurement data, the frequency correlation coefficient initial value is 0.2-0.5, the polarizability initial value is that the time constant initial value is the inverse of the frequency of relative phase polarographic maximum correspondence with the maximum frequency dispersion rate of spectral amplitude definition;

3) use MCMC method (randomness algorithm): the initial value that the model optimum solution of being calculated employing step 2) is extracted as the MCMC method, when being provided with, MCMC method parameter provides parameter value priori scope to be measured according to the actual conditions in image data or the application process, and in computing machine, generate the multiple parameter model of Bayes at random by the foundation of priori scope, produce and analyze the Metropolis-Hasting Markov chain with the MCMC method by the Gaussian distribution method then, thereby try to achieve the overall edge posterior probability density function of each unknown parameter in effective range; According in the step 1) the given standard accepted check generation model and calculate the extraction precision;

Described parameter value priori scope to be measured is: resistivity priori scope is a complex resistivity spectral amplitude variation range, frequency correlation coefficient priori value range is 0.001-0.99, polarizability priori scope is 0.001-0.99, and time constant priori scope is 0.001-10000s;

4) testing model: the gained precision does not reach the regulation requirement if extract as a result, and the overall probability density function that extract with MCMC method in step 3) this moment asks the parameter initial value still to carry out step 3) as waiting, so repeats to meet the demands up to extracting precision; If extracting precision reaches requirement and then uses step 2) in the overall probability density function that extracts with the MCMC method of this moment of least square method be initial value, thereby extract the globally optimal solution that obtains on the practical significance again;

5) result to be extracted stable after, carry out match to data and the result analyzed extracting curve and institute at last.

Images