CN101706587B - 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 method for distilling 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\ω}{\mathrm{\τ}}_1\right)}^{c_1}}]-m_2[1-\frac{1}{1+{\left(\mathrm{i\ω}{\mathrm{\τ}}_2\right)}^{c_2}}]+m_3[1-\frac{1}{1+{\left(\mathrm{i\ω}{\mathrm{\τ}}_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 it and exist, though can address this problem to a certain extent through the model of above-mentioned Cole-Cole frequency dispersion sum; 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, and the degree of accuracy that the data demand that it is also surveyed instrument is higher.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, and accuracy of instrument is had relatively high expectations, and 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 objective of the invention is purpose through following technical scheme realization:

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 representes frequency ratio, is constant, k＞1.If adopt 2 ⁿSeries pseudorandom multi-frequency signal is as excitation field source, and 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 EM coupling, can suppress inductive coupling; The observation of relative phase spectrum is measured than frequency dispersion rate has 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 through 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 through the extraction to single Cole-Cole model relative phase spectrum and spectral amplitude, single complex resistivity Cole-Cole expression formula is following:

$\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="H2009102270095E00031.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="H2009102270095E00031.png"\; state="\{\{state\}\}">c</figure-callout>}}{2}+{(1-m)}^2{\left(\mathrm{\&omega;\&tau;}\right)}^{2\mathrm{<figure-callout\; id="1"\; label="c"\; filenames="H2009102270095E00011.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="H2009102270095E00031.png"\; state="\{\{state\}\}">c</figure-callout>}}{2}+{\left(\mathrm{\&omega;\&tau;}\right)}^{2c}}\right]}^{1/2}$

Its phase spectrum

is:

$=\mathrm{arctg}\frac{-\mathrm{mI}}{(1-m)(R^2+I^2)+\mathrm{mR}}---\left(4\right)$

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

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

$-\frac{1}{k-1}\frac{-{\mathrm{\&rho;}}_{0}m}{{\left[A\left(\mathrm{k\&omega;}\right)\right]}^2}[\frac{({\mathrm{<figure-callout\; id="2"\; label="I"\; filenames="H2009102270095E00031.png"\; state="\{\{state\}\}">I</figure-callout>}}_2^2-{\mathrm{<figure-callout\; id="2"\; label="R"\; filenames="H2009102270095E00031.png"\; state="\{\{state\}\}">R</figure-callout>}}_2^2)\mathrm{Im}\left(\mathrm{k\&omega;}\right)-2R_2{\mathrm{<figure-callout\; id="2"\; label="I"\; filenames="H2009102270095E00031.png"\; state="\{\{state\}\}">I</figure-callout>}}_2\&CenterDot;\mathrm{Re}\left(\mathrm{k\&omega;}\right)}{{({\mathrm{<figure-callout\; id="2"\; label="R"\; filenames="H2009102270095E00031.png"\; state="\{\{state\}\}">R</figure-callout>}}_2^2+{\mathrm{<figure-callout\; id="2"\; label="I"\; filenames="H2009102270095E00031.png"\; state="\{\{state\}\}">I</figure-callout>}}_2^2)}^2}\&CenterDot;\frac{\&PartialD;{\mathrm{<figure-callout\; id="2"\; label="R"\; filenames="H2009102270095E00031.png"\; state="\{\{state\}\}">R</figure-callout>}}_2}{\&PartialD;c}-\frac{2R_2{\mathrm{<figure-callout\; id="2"\; label="I"\; filenames="H2009102270095E00031.png"\; state="\{\{state\}\}">I</figure-callout>}}_2\&CenterDot;\mathrm{Im}\left(\mathrm{k\&omega;}\right)+({\mathrm{<figure-callout\; id="2"\; label="I"\; filenames="H2009102270095E00031.png"\; state="\{\{state\}\}">I</figure-callout>}}_2^2-{\mathrm{<figure-callout\; id="2"\; label="R"\; filenames="H2009102270095E00031.png"\; state="\{\{state\}\}">R</figure-callout>}}_2^2)\&CenterDot;\mathrm{Im}\left(\mathrm{k\&omega;}\right)}{{({\mathrm{<figure-callout\; id="2"\; label="R"\; filenames="H2009102270095E00031.png"\; state="\{\{state\}\}">R</figure-callout>}}_2^2+{\mathrm{<figure-callout\; id="2"\; label="I"\; filenames="H2009102270095E00031.png"\; state="\{\{state\}\}">I</figure-callout>}}_2^2)}^2}\&CenterDot;\frac{\&PartialD;{\mathrm{<figure-callout\; id="2"\; label="I"\; filenames="H2009102270095E00031.png"\; state="\{\{state\}\}">I</figure-callout>}}_2}{\&PartialD;c}]---\left(6\right)$

$-\frac{1}{k-1}\frac{-{\mathrm{\&rho;}}_{0}m}{{\left[A\left(\mathrm{k\&omega;}\right)\right]}^2}[\frac{({\mathrm{<figure-callout\; id="2"\; label="I"\; filenames="H2009102270095E00031.png"\; state="\{\{state\}\}">I</figure-callout>}}_2^2-{\mathrm{<figure-callout\; id="2"\; label="R"\; filenames="H2009102270095E00031.png"\; state="\{\{state\}\}">R</figure-callout>}}_2^2)\mathrm{Im}\left(\mathrm{k\&omega;}\right)-2R_2{\mathrm{<figure-callout\; id="2"\; label="I"\; filenames="H2009102270095E00031.png"\; state="\{\{state\}\}">I</figure-callout>}}_2\&CenterDot;\mathrm{Re}\left(\mathrm{k\&omega;}\right)}{{({\mathrm{<figure-callout\; id="2"\; label="R"\; filenames="H2009102270095E00031.png"\; state="\{\{state\}\}">R</figure-callout>}}_2^2+{\mathrm{<figure-callout\; id="2"\; label="I"\; filenames="H2009102270095E00031.png"\; state="\{\{state\}\}">I</figure-callout>}}_2^2)}^2}\&CenterDot;\frac{\&PartialD;{\mathrm{<figure-callout\; id="2"\; label="R"\; filenames="H2009102270095E00031.png"\; state="\{\{state\}\}">R</figure-callout>}}_2}{\&PartialD;\mathrm{\&tau;}}-\frac{2R_2{\mathrm{<figure-callout\; id="2"\; label="I"\; filenames="H2009102270095E00031.png"\; state="\{\{state\}\}">I</figure-callout>}}_2\&CenterDot;\mathrm{Im}\left(\mathrm{k\&omega;}\right)+({\mathrm{<figure-callout\; id="2"\; label="I"\; filenames="H2009102270095E00031.png"\; state="\{\{state\}\}">I</figure-callout>}}_2^2-{\mathrm{<figure-callout\; id="2"\; label="R"\; filenames="H2009102270095E00031.png"\; state="\{\{state\}\}">R</figure-callout>}}_2^2)\&CenterDot;\mathrm{Re}\left(\mathrm{k\&omega;}\right)}{{({\mathrm{<figure-callout\; id="2"\; label="R"\; filenames="H2009102270095E00031.png"\; state="\{\{state\}\}">R</figure-callout>}}_2^2+{\mathrm{<figure-callout\; id="2"\; label="I"\; filenames="H2009102270095E00031.png"\; state="\{\{state\}\}">I</figure-callout>}}_2^2)}^2}\&CenterDot;\frac{\&PartialD;I_2}{\&PartialD;c}]---\left(7\right)$

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

$\left.\begin{array}{c}{\mathrm{<figure-callout\; id="1"\; label="R"\; filenames="H2009102270095E00011.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="H2009102270095E00011.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="H2009102270095E00031.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="H2009102270095E00031.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="H2009102270095E00031.png"\; state="\{\{state\}\}">c</figure-callout>}}{2}\end{array}\right\}$

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

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

${\mathrm{<figure-callout\; id="2"\; label="I"\; filenames="H2009102270095E00031.png"\; state="\{\{state\}\}">I</figure-callout>}}_2^2-{\mathrm{<figure-callout\; id="2"\; label="R"\; filenames="H2009102270095E00031.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="H2009102270095E00031.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="H2009102270095E00031.png"\; state="\{\{state\}\}">c</figure-callout>}}{2}$

$\frac{\&PartialD;{\mathrm{<figure-callout\; id="1"\; label="R"\; filenames="H2009102270095E00011.png"\; state="\{\{state\}\}">R</figure-callout>}}_1}{\&PartialD;c}={\left(\mathrm{\&omega;\&tau;}\right)}^c[\mathrm{in}\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="H2009102270095E00031.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="H2009102270095E00031.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="H2009102270095E00031.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="H2009102270095E00031.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="H2009102270095E00011.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="H2009102270095E00031.png"\; state="\{\{state\}\}">c</figure-callout>}}{2}$ $\frac{\&PartialD;{\mathrm{<figure-callout\; id="2"\; label="R"\; filenames="H2009102270095E00031.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="H2009102270095E00031.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="H2009102270095E00031.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="H2009102270095E00031.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 following to draw its extraction fundamental formular according to formula (3):

$\frac{\&PartialD;A}{\&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="H2009102270095E00031.png"\; state="\{\{state\}\}">c</figure-callout>}}{2}+{(1-m)}^2{\left(\mathrm{\&omega;\&tau;}\right)}^{2\mathrm{<figure-callout\; id="1"\; label="c"\; filenames="H2009102270095E00011.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="H2009102270095E00031.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 method for distilling that is appropriate to the relative phase spectrum.The technical scheme that the present invention adopted is to combine random statistical algorithm MCMC method (the special calot's method of Marko Fumeng) to unite with the least square method in the deterministic algorithm 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 through polynary Gaussian distribution method.

Whole technical proposal may further comprise the steps:

1) go out with the Cole-Cole Model Calculation of giving that relative phase is composed and measured value that spectral amplitude theoretical value or directly use has recorded in the instrument as the data of model extraction; Selected Cole-Cole model parameter to be asked waits to ask the theory relation between model parameter and relative phase spectrum and spectral amplitude observation data according to extraction model formula (5), (6), (7), (9) foundation; 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): on original relative phase spectrum and amplitude spectrum model based, setting up the least square model is least square and model; And set by step 1) given standard provides relevant parameter setting in; Adopt relative phase spectrum and amplitude spectrum data given in the step 1); Appoint based on 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 said parameter initial value is: the resistivity initial value is the corresponding resistivity of spectral amplitude low-limit frequency in the measurement data; The frequency correlation coefficient initial value is 0.2-0.5; The polarizability initial value is the maximum frequency dispersion rate with the spectral amplitude definition, and the time constant initial value is the inverse of the corresponding frequency of relative phase polarographic maximum.

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;, MCMC method parameter provides parameter value priori scope to be measured when being provided with 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 through 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;

Said 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 with institute to data and the result analyzed extracting curve at last.

The present invention extracts 4 parameters of Cole-Cole model through 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, thus 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 basically.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 through MCMC method and least square method associating method for distilling that not choosing of initial value can impact and obtain a large amount of extraction uncertainty information to extracting the result, also retained least square method extraction rate advantage faster simultaneously, proved by our experiment; This associating method for distilling 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

Below in conjunction with accompanying drawing and embodiment the present invention is further specified.

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

1. the theoretical value and the calculated value behind the relative noise of adding of calculating them with the formula of existing Cole-Cole model relative phase spectrum and spectral amplitude according to 20 given class frequency values are organized data 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 establish an initial value in the parameter value scope, the least square method algorithm that utilizes the front to provide is obtained the one of which local minimum and is begun the extraction initial value as the MCMC method; The extraction process step of least square method (Gauss-Newton method) is following:

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=1,2 ..., M; J=1,2 ..., N)

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

3) group of solving an equation: $A_k^T{A}_k{b}^{\left(k\right)}=-A_k^T{f}^k,$ Get Gauss-newton's direction p ^(k)

4) from b ^(k)Set out, along direction p ^(k)Do linear search, confirm 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 like ρ ₀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 parameter has been set, 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 following that the MCMC method is extracted idiographic flow:

(1) try to achieve the likelihood equation that relative phase is composed through Gaussian distribution:

μ 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 let 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) let 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 needing to extract extracts as the initial value in the 2nd step once again 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 to stop to extract when perhaps its standard deviation that calculates is within error range after the extraction result is stable.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 practical implementation step is following:

1. spectral amplitude and relative phase spectral theory calculated value are extracted, can inspection MCMC method reach our required extraction accuracy requirement; So need obtain another parameter ρ through 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 through 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; Carry out this experiment in step and be in order to detect relative phase spectrum and spectral amplitude in the extraction precision that has noise jamming under with the MCMC method, theoretical value is still continued to use the 1st value in going on foot; 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.See that from extracting curve and result except time constant, other parameter extraction all compares near actual value, 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 EM coupling in theoretical value is calculated; Extract with single Cole-Cole model relative phase spectrum then, through extracting the compacting electromagnetic coupling effect effect that the result confirms that relative phase is composed.The EM coupling 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, the extraction error is 1.0889e-004.Fig. 8 is τ=10 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, the extraction error is 6.5350e-005.Fig. 9 is τ=10 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, the extraction error is 2.3529e-004.From the extraction result of 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 is shown in figure 10, and it extracts the result and is: τ=7.29688e+010; C=0.151606; M=0.916147, extracting error is 0.000014, and owing to locate unusually at low-resistivity, the time constant value has been difficult to accurate Calculation, and therefore above-mentioned data can only reflect that low-resistivity is that copper coin is unusual unusually, and c and m approach actual value.

Above explanation mode only is used to explain the present invention; And be not the restriction to protection domain of the present 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 method for distilling 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 through 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: combine the least square method in random statistical algorithm MCMC method and the deterministic algorithm, from relative phase spectrum and spectral amplitude, extract the Cole-Cole model parameter; 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 through polynary Gaussian distribution method;

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

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

In said 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\}$

$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}$

$I_1^2-R_1^2={\left(\mathrm{\&omega;\&tau;}\right)}^{2c}({\sin }^2\frac{\mathrm{\&pi;c}}{2}-{\cos }^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}={\mathrm{k\&omega;\&tau;}}^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 on original relative phase spectrum and amplitude spectrum model based, to set up the least square model is least square and model; And set by step 1) given standard provides relevant parameter setting in; Adopt relative phase spectrum and amplitude spectrum data given in the step 1); Appoint based on the scope of waiting to ask parameter 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 said parameter initial value is: the resistivity initial value is the corresponding resistivity of spectral amplitude low-limit frequency in the measurement data; The frequency correlation coefficient initial value is 0.2-0.5; The polarizability initial value is the maximum frequency dispersion rate with the spectral amplitude definition, and the time constant initial value is the inverse of the corresponding frequency of relative phase polarographic maximum;

3) use randomness algorithm MCMC method: the initial value that the model optimum solution of being calculated employing step 2) is extracted as the MCMC method;, MCMC method parameter provides parameter value priori scope to be measured when being provided with 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 through 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;

Said 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 with institute to data and the result analyzed extracting curve at last.

Images