CN108169801A - High-resolution earth resistivity fast imaging method - Google Patents
High-resolution earth resistivity fast imaging method Download PDFInfo
- Publication number
- CN108169801A CN108169801A CN201810038410.3A CN201810038410A CN108169801A CN 108169801 A CN108169801 A CN 108169801A CN 201810038410 A CN201810038410 A CN 201810038410A CN 108169801 A CN108169801 A CN 108169801A
- Authority
- CN
- China
- Prior art keywords
- resistivity
- depth
- gds
- curves
- curve
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V3/00—Electric or magnetic prospecting or detecting; Measuring magnetic field characteristics of the earth, e.g. declination, deviation
- G01V3/38—Processing data, e.g. for analysis, for interpretation, for correction
Landscapes
- Life Sciences & Earth Sciences (AREA)
- Engineering & Computer Science (AREA)
- Environmental & Geological Engineering (AREA)
- Geology (AREA)
- Remote Sensing (AREA)
- Physics & Mathematics (AREA)
- General Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Geophysics (AREA)
- Geophysics And Detection Of Objects (AREA)
- Measurement Of Resistance Or Impedance (AREA)
Abstract
The present invention relates to a kind of high-resolution earth resistivity fast imaging methods, inverting and high-resolution earth resistivity fast imaging including high-resolution earth resistivity image reconstruction, the former is on the basis of enlightening inversion method is routinely helped, it differentiates to the resistivity curve that inverting obtains, GDS ~ Z figures are drawn, realize differential imaging；The latter is directed to three layers of earth-electricity model, on the basis of the Da Zhanuoke methods of inversion, respectively four kinds of difference Da Zhanuoke inversional curves are corrected with the fitting of coefficient, the correction coefficient that first layer and the second layer are fitted by nonlinear least square fitting function 1squcurvefit draws GDS ~ Z figures, realizes differential imaging.The present invention can compare the position for being intuitive to see lithology interface, using GDT high-resolution survey meter detection datas as research object, by GDS differential imaging technologies, improve the resolution capability to thin layer, the precision and reliability of inverting are obtained for raising.
Description
Technical field
The present invention relates to Engineering geophysical exploration fields, and in particular to a kind of high-resolution earth resistivity fast imaging side
Method.
Background technology
Resistivity imaging technology is always the advanced topic that geophysics circle is paid close attention to and most complicated, most
The new and high technology of challenge, it is special that resistivity imaging technology has in a series of problems, such as solving resource, mineral products, environment and engineering
Meaning.1987, Shima was put forward for the first time " resistivity tomography " word, and the method for proposing inversion interpretation, this
Afterwards, many people, particularly Japan and the United States scholar have carried out research from theoretical, experiment to this problem to application from different perspectives.Nearly ten
With the development of array cloth pole mode data collecting system over year, it is already possible to collect potential data intensive enough, have
The condition of imaging method, domestic scholars started in the last century 90's this respect research (Wang Yan equalitys, 1994；Week
Soldier etc., 1995；Wang Xingtai etc., 1995；Bottom high official position etc., 1997；Yan Zhongqiong etc., 1996；Feng Rui etc., 1997；Li Zhi is bright, and 2003)
A small amount of tentative application study is completed, as a kind of new technology, resistivity tomography reaches far away the stage of ripeness, is being permitted
It is all left to be desired on more sport technique segments.There are the problem of：1. majority resistivity method is all finding better observed pattern at present
To adapt to high-precision, high density, high-resolution imaging.2. for intricately electric structure, there are many model parameter, global inversion algorithm
It is being currently unrealistic for applying.3. there are low-resistance to " the intrusion effect " of high resistant in some imaging methods, with depth
Increasing resolution ratio reduces and occurs with abnormal.
In recent years, people did a lot on the method and technology of resistivity imaging, propose some observation procedures and
Algorithm, development trend are mainly reflected in the following aspects：
1st, the improvement of observed pattern, to adapt to the needs of imaging, GDS methods are exactly in recent years, to be designed by the Ministry of Railways first
A kind of observed pattern for " three high " (high-precision, high density, high-resolution) that Xi’an Branch of institute proposes is laid a good foundation for imaging.
2nd, the improvement of imaging method
1. Zohdy method
The method is that a kind of of iteration and fitting method represents method, actually a kind of Least-squares minimization method, the novelty of the method
Part is that its initial model and how to be exported that the electrical image that this method obtains is valuable in geology.
Especially when there is quality data, the superiority of method can be more showed, by the GDS methods with " three high " technology and assistant enlightening inverting
Method is combined, and the advantages of being a wise selection, making its method more highlightedly embodies.
2. Electric current trace method
The method is ray tracing technique and the Tracing current ray current potential resistivity imaging side that proposes when being walked in analogy seismology
Method, not only related with the resistivity of the unit by the potential difference of each unit in glucose current equation line tracing process, Er Qiehe
Other resistivity distributions are related.Therefore, there has been proposed a kind of improved two spots ray tracing technologies, the wherein meter of potential difference
The variation in view of current density is calculated, is divided into two levels, first, square being inversely proportional due to geometrical attenuation and distance；Second is that by
It is nonlinear in the law of refraction of current line, need to be corrected by refraction theorem.Therefore, although the method, which is one, extremely hair
Open up future method, but there are the problem of it is quite a few, be a theory and apply the upper method all needed to be further improved.
3. equipotential line tracing method
The imaging method that the method has been used for reference current potential imaging (APT) technology medically applied and proposed, basic principle are：
The current potential measured is tracked from measurement point along equipotential line to underground, the resistivity value of each underground unit that modification equipotential line is passed through,
A simple projection model is can obtain, but there are one condition, exactly it has to be assumed that electric current is limited to these etc.
It is flowed in bit line narrowband, it is clear that this is a simple and quick algorithm, but also has shortcoming, such as low-resistance " invading to high resistant
Enter effect ", increasing to differentiate with depth reduces and occurs with abnormal.
By above to the summary of resistivity tomography method as can be seen that resistivity CT is studied at nearly twenties years
Achieve many important achievements in development, but due to the specific condition of geophysics, in data critical shortage, the earth
Portion is complicated etc. so that theory can solve the problems, such as at present and objective reality situation also has larger gap.Moreover, to the greatest extent
Managing various methods respectively has feature, the cross-section image of subsurface resistivity can be obtained, but resolution ratio is not very high, with actual production
Using also a certain distance；More to the research of algorithm simultaneously, the research to method is not too deep.
Invention content
The object of the present invention is to provide a kind of high-resolution earth resistivity fast imaging methods, i.e. GDS differential imaging methods, are used for
It solves the problems, such as that conventional earth resistivity imaging is inaccurate and precision is inadequate, is effectively improved imaging precision and reliability, explains
As a result it is more intuitive..
The technical solution adopted in the present invention is：
High-resolution earth resistivity fast imaging method, it is characterised in that：
Include the following steps：
(1) inverting of high-resolution earth resistivity image reconstruction on the basis of routinely enlightening inversion method is helped, obtains inverting
Resistivity curve differentiate, draw GDS~Z figure, realize differential imaging；
(2) high-resolution earth resistivity fast imaging, for three layers of earth-electricity model, on the basis of the Da Zhanuoke methods of inversion,
Respectively four kinds of difference Da Zhanuoke inversional curves are corrected with the fitting of coefficient, passes through nonlinear least square fitting function
1squcurvefit fits first layer and the correction coefficient of the second layer draws GDS~Z figures, realizes differential imaging.
The step of inverting of high-resolution earth resistivity image reconstruction, includes：
Establish initial model；
Determine the depth of each layer；
Using Forward Formula, resistivity transforming function transformation function T is calculatedL(k)；
Build minimum object function δ；
It obtains meeting the minimum optimum stepsize vector Δ P of δ, then to PiIt is corrected, acquires Pi+1Point (Pi+1=Pi+Δ
P)；
It iterates, Approach by inchmeal makes TL(k,Pi) really close to T1(k), more it is accurate to Δ P so as to acquire；
Inversion result is obtained, draws ρ (z)~z curves；
Ask for GDS values；
Draw GDS~z curves namelyCurve；
Data are carried out with multiple structure using Smoothing Technique and Cubic Spline Interpolation is slided；Then, to reconstructing later number
According to high-density sampling is carried out, differential imaging is finally realized.
The step of high-resolution earth resistivity fast imaging, includes：
ρs(r)~r is seen as Da Zhanuoke curves
It calculates
Calculate depth z and corresponding resistivity；
Paint ρ (z)~z curves；
Build the relational expression of adaptation coefficient；
Export adaptive correction coefficient；
ρ (z)~z curves are corrected by correction coefficient, regenerate ρ (z)~z curves；
Ask for GDS values；
Draw GDS~z curves namelyCurve；
Data are carried out with multiple structure using Smoothing Technique and Cubic Spline Interpolation is slided；Then, to reconstructing later number
According to high-density sampling is carried out, GDS differential imagings are finally realized.
The specific steps of the inverting of high-resolution earth resistivity image reconstruction include：
Step 1：Establish initial model：
Assuming that stratum the number of plies and sounding curve points as many, every layer of resistivity is respective points on sounding curve
Apparent resistivity, every layer of depth is equal to the electrode spacings of respective points on sounding curve multiplied by with a depth shift factor, this
Sample has just obtained an initial model (hi,ρi, i=1,2 ... is n)；
Step 2：Determine the depth of each layer：
Since in inversion algorithm, the depth initial value of each layer is to be replaced with electrode spacing multiplied by with a depth shift factor
, i.e.,To accelerate inverting convergence rate in inverting, following iteration is used to the inversion result of depth
Form：
H(i+1)=(H(i)+ΔH(i))·K (1)
Wherein K is the depth shift factor；The selection method of K is attributed to：
1. setting initial the depth shift factor as K=1, the depth value of each layer of initial model is equal to corresponding electrode spacing value；
2. calculate the theoretical apparent resistivity value ρ of initial modelC；
3. calculate theoretical apparent resistivity and the root-mean-square error R of measured visual resistivity depth measurement data：
Wherein：I is iterations, and j is jth layer or j-th of electrode spacing, ρ0(j) it is that actual measurement on j-th of electrode spacing regards electricity
Values of resistivity, ρci(j) the theoretical apparent resistivity value to be calculated on j-th of electrode spacing of iteration j；
If 4. R0>R1, then the depth shift factor, i.e. K=0.9K are reduced, 2. 3. repetition walks, and compare root-mean-square error, directly
To Ri>Ri-1When, the depth shift factor K at this moment is optimum value；
Step 3：Using Forward Formula, resistivity transforming function transformation function T is calculatedL(k)；
Step 4：Build minimum object function δ；
The T that forward modeling is obtainedL(k) the transforming function transformation function T with surveying1(k) compare, the variance of the two is obtained；
M is T in formula1(k) hits, it is equal to ρs(r) hits, and
I is iterations, i=0,1,2 ...；Model as i=0 is initial model；P is vector；N is the layer of model
Number；
If by TL(k,Pi) in the initial value P of ith iterationiNearby it is unfolded by Tailor progression, secondary above height will be omitted
Rank, then
So as to obtain：
Number of parameters of the N=2n-1 for layer, Δ P in formulaiIt is j-th of parameter PjCorrection amount or correct step-length；
Step 5：According to the extreme value theory of the function of many variables, formula (5) is converted, can be obtained：
Step 6：Above-mentioned equation is solved, just obtains meeting the minimum optimum stepsize vector Δ P of δ；Then to PiIt is corrected, asks
Obtain Pi+1Point (Pi+1=Pi+ΔP)；
Step 7：It repeats above-mentioned step to gather, iterate, Approach by inchmeal makes TL(k,Pi) really close to T1(k), so as to ask
It obtains and is more accurate to Δ P；
In inverting iteration, in order to improve iteration speed, each iteration with measured visual resistivity and is calculated regarding electricity as the following formula
The ratio between resistance rate is coefficient, adjusts the resistivity of each layer：
Specific practice is：Every time in iterative calculation, root-mean-square error R is judgediWhether it is less than ε, takes ε=5%, if full
Foot requirement, at this moment model parameter resistivity, depth be inverting as a result, otherwise by above formula adjust resistivity, then sentence again
It is disconnected whether to have R2<ε repeats above step and gathers, until equal root error is met the requirements；
Step 8：Inversion result is obtained, draws ρ (z)~z curves；
Step 9：Ask for GDS values；GDS is asked for according to formula (9)：
Wherein, ziFor i-th of depth value, zi-1For (i-1)-th depth value, resistivity values of the ρ (i) for i-th of depth, ρ
(i-1) resistivity value for (i-1)-th depth；
Step 10：Draw GDS~z curves namelyCurve；
By rightThe identification of curvilinear characteristic can more intuitively reflect the variation tendency of electrical parameter；
Step 11：Due to being limited by electric sounding data, measured line smoothing degree is inadequate, is unfavorable at differential imaging
Reason；Therefore, in order to preferably hold the feature of curve consecutive variations and accurate derivation；Using slip Smoothing Technique and cubic spline
Function interpolation carries out data multiple structure；Then, high-density sampling is carried out to reconstructing later data, finally realizes differential imaging.
The specific steps of high-resolution earth resistivity fast imaging include：
Step 1：ρs(r)~r is seen as Da Zhanuoke curves
Step 2：It calculates
Step 3：It calculates depth z and corresponding resistivity, specific formula for calculation is as follows：
Step 4：Paint ρ (z)~z curves；
Step 5：Build the relational expression of adaptation coefficient；
For three layers of earth-electricity model, if H1And H2It is the first bed boundary of earth-electricity model and the depth of the second bed boundary respectively,
F1And F2It is that corresponding earth-electricity model carries out the first bed boundary and the second bed boundary that are obtained after Da Zhanuoke curve method invertings respectively
Depth；It analyzes through a large number of experiments, finds H1/F2And F2/F1Between, H2/F1And F1/F2Between there are certain relationship, take
X=F1/F2；For by the obtained initial results of Da Zhanuoke curve method invertings, the H of four type electrical sounding curves1/F2With
F2/F1Between functional relation be all nonlinear, it is possible to nonlinear least square fitting function lsqucurvefit come
It is fitted both the above functional relation；
Step 6：Adaptive correction coefficient exports；
For the correction coefficient of the first bed boundary,It is corrected for the second bed boundary
Coefficient, a2=x*f (x).Through with nonlinear least square fitting method to H1/F2With F2/F1And H2/F1With F1/F2Distribution function
After fitting, the fitting function f (x) and interface correction coefficient a of four kinds of earth-electricity models can be respectively obtained；
For H-type geoelectric cross section：
First bed boundary fitting function：
First bed boundary correction coefficient：
Second bed boundary fitting function：
f2(x)=31.0405*exp (- 43.6255*x)+4.7257*exp (- 5.5429*x)
(13)
Second bed boundary correction coefficient：
Step 7：The ρ (z) in step 4~z curves are corrected by the correction coefficient of step 6, regenerate ρ (z)
~z curves；
Step 8：GDS values are asked for, GDS is asked for according to formula (9)；
Step 9：Draw GDS~z curves namelyCurve；
By rightThe identification of curvilinear characteristic can more intuitively reflect the variation tendency of electrical parameter；
Step 10：Due to being limited by electric sounding data, measured line smoothing degree is inadequate, is unfavorable at differential imaging
Reason；Therefore, in order to preferably hold the feature of curve consecutive variations and accurate derivation；Using slip Smoothing Technique and cubic spline
Function interpolation carries out data multiple structure；Then, carry out high-density sampling to reconstructing later data, finally realize GDS differential into
Picture.
In step 1-4：
1. on first branch asymptote, ρs(r)=ρ1, a=0, therefore have ρ (z)=ρ1；
2. in ρnFor on limited tail branch asymptote, ρs(r)=ρn, a=0, so ρs(r)=ρn；
3. work as ρnDuring → ∞, the slope of tail branch asymptote is 45 ° of straight lines, at this time a=1, therefore ρ (z)=∞；
4. work as ρnWhen → 0, the slope a of tail branch asymptote<- 1 will lose meaning；Work as a during practical calculating<When -1, a should be taken
=-1.
The present invention has the following advantages：
(1) to the resolution ratio higher of thin layer, depth divides more accurate (such as Fig. 3,4).
(2) by a large amount of Numerical modelling, three layers of earth-electricity model are given to different types of Da Zhanuoke curves
Correction coefficient.Result of calculation show it is corrected after as a result, having large increase on Explanation Accuracy, in normal condition
Under, the error of inversion result and earth-electricity model is within ± 5%, and the characteristics of method is that calculating speed is fast, as a result of adaptive
Bearing calibration is answered, human factor influences small.
(3) GDS inversion imagings method has the characteristics of prominent weak anomaly, high resolution compared with conventional electric sounding inverting.
(4) successful application of GDS differential imagings so that explanation results are more intuitive, with the increasing of imaging data density,
GDS imaged images will will appear.
Description of the drawings
Fig. 1 is the Forward Modeling and Inversion of K-type and its first derivative figure.(ρ1=10 Ω .m, ρ2=60 Ω .m, ρ3=5 Ω .m, h1=
19m,h 2=209m；- --- is bent for derivative；--- --- is ρ (H)-H curves；For forward modeling curve).
Fig. 2 is the Forward Modeling and Inversion of Q types and its first derivative figure.(ρ1=100 Ω .m, ρ2=10 Ω .m, ρ3=100 Ω .m, h1=15m, h2=120m.)
Fig. 3 is the good model for leading thin layer, forward modeling curve and GDS curves.(ρ1=100 Ω .m, ρ2=20 Ω .m, ρ3=100
Ω.m,h 1=100m, h2=3m.)
Fig. 4 is model, forward modeling curve and the GDS curves of high resistant thin layer.(ρ1=50 Ω .m, ρ2=200 Ω .m, ρ3=50
Ω.m,h1=100m, h2=3m.)
Fig. 5 rebuilds inverting and forward modeling Comparative result for Q type sections resistivity image.(dotted line is forward modeling curve, and solid line is anti-
Drill curve, ρ1=100 Ω .m, ρ2=40 Ω .m, ρ3=10 Ω .m, h1=20m, h2=100m)
Fig. 6 is G type section GDS imaging results.(ρ1=10 Ω .m, ρ2=100 Ω .m, h=20m)
Fig. 7 is H-type section GDS imaging results.(ρ1=100 Ω .m, ρ2=10 Ω .m, ρ3=100 Ω .m, h1=20m, h2
=100m)
Fig. 8 is K-type section GDS imaging results.(ρ1=10 Ω .m, ρ2=60 Ω .m, ρ3=20 Ω .m, h1=20m, h2=
100m)
Fig. 9 is AA type section GDS imaging results.(ρ1=10 Ω .m, ρ2=20 Ω .m, ρ3=60 Ω .m, ρ4=80 Ω .m,
h1=20m, h2=60m, h3=100m)
Figure 10 is H-type iunction for curve distribution map.(ρ1=500 Ω .m, ρ2=50 Ω .m, ρ3=450 Ω .m, h1=
20m,h2=200m)
For H-type the electric disconnected GDS imaging results of Figure 11.(ρ1=500 Ω .m, ρ2=50 Ω .m, ρ3=450 Ω .m, h1=
20m,h2=200m)
For K-type the electric disconnected GDS imaging results of Figure 12.(ρ1=5 Ω .m, ρ2=50 Ω .m, ρ3=5 Ω .m, h1=12m, h2=
212m)
For A types the electric disconnected GDS imaging results of Figure 13.(ρ1=5 Ω .m, ρ2=50 Ω .m, ρ3=400 Ω .m, h1=17m, h2
=160m)
For Q types the electric disconnected GDS imaging results of Figure 14.(ρ1=200 Ω .m, ρ2=30 Ω .m, ρ3=4 Ω .m, h1=8m, h2
=180m)
Specific embodiment
The present invention will be described in detail With reference to embodiment.
The present invention relates to a kind of high-resolution earth resistivity fast imaging method, including high-resolution earth resistivity image reconstruction
Inverting and high-resolution earth resistivity fast imaging two parts content.
GDS methods with " three high " technology are combined by first part with traditional assistant enlightening method of inversion.Using GDS differential into
As technology, it can more intuitively reflect the change shape of electrical parameter, the extreme value pair of the first derivative curve of ρ (z)~z curves
Answer the boundary position (as shown in Figure 1) of lithologic interface.Inversion result has the characteristics that reliable, with high accuracy, eliminates low-resistance and " invades
Enter effect " influence.
GDS methods with " three high " technology are combined by second part with the Da Zhanuoke methods of inversion.And in Zha Nuoke invertings
On the basis of carried out a large amount of improvement, by a large amount of Numerical modelling, different types of Da Zhanuoke curves are provided
The correction coefficient of three layers of earth-electricity model.Using GDS differential imaging technologies, it can more intuitively reflect the variation of electrical parameter
Form.Effectively increase the precision and reliability of inverting.
High-resolution Electric Sounding Technology in Survey of Wellpit (GDS) is realized as a result of GDT high-resolution survey meter and its distinctive observation technology
The high-density sampling of the quick multi-fold tolling measurement of equal difference pole span, instrumental sensitivity is high, and stability is good, can be in low background
Lower extraction weak anomaly.Research is unfolded in GDS differential imaging methods on the basis of high-resolution Electric Sounding Technology in Survey of Wellpit (GDS), with
Conventional electric sounding inverting is compared, and is had the characteristics of prominent weak anomaly, high resolution, is overcome conventional inversion method for height group
The problem of thin layer depth determines inaccuracy (such as Fig. 3,4) substantially increases the precision and reliability of inversion result.
Based on conditions above, to amplify, protruding local anomaly, low noise Anomalies of Backgrounds is suppressed, to the important ginseng of this method
Number GDS is defined as follows：
Wherein, ziFor i-th of depth value, zi-1For (i-1)-th depth value, resistivity values of the ρ (i) for i-th of depth, ρ
(i-1) resistivity value for (i-1)-th depth.
In order to ensure that GDS has higher resolution ratio, Δ z is generally required<<H (H is destination layer buried depth), generally takes
Specific implementation step of the present invention is as follows：
First part：The inversion method of high-resolution earth resistivity image reconstruction
This method is a kind of Least-squares minimization method of inversion proposed for GDS methods.Step includes：
Step 1：Establish initial model.
Assuming that stratum the number of plies and sounding curve points as many, every layer of resistivity is respective points on sounding curve
Apparent resistivity, every layer of depth is equal to the electrode spacing of respective points on sounding curve multiplied by with a depth shift factor.This
Sample has just obtained an initial model (hi,ρi, i=1,2 ... is n).
Step 2：Determine the depth of each layer.
Since in inversion algorithm, the depth initial value of each layer is to be replaced with electrode spacing multiplied by with a depth shift factor
, i.e.,To accelerate inverting convergence rate in inverting, following iteration is used to the inversion result of depth
Form：
H(i+1)=(H(i)+ΔH(i))·K (2)
Wherein K is the depth shift factor.And the selection method of K can be attributed to：
1. setting initial the depth shift factor as K=1, the depth value of each layer of initial model is equal to corresponding electrode spacing value.
2. calculate the theoretical apparent resistivity value ρ of initial modelC。
3. calculate theoretical apparent resistivity and the root-mean-square error R of measured visual resistivity depth measurement data：
Wherein：I is iterations, and j is jth layer or j-th of electrode spacing, ρ0(j) it is that actual measurement on j-th of electrode spacing regards electricity
Values of resistivity, ρci(j) the theoretical apparent resistivity value to be calculated on j-th of electrode spacing of iteration j.
If 4. R0>R1, then the depth shift factor, i.e. K=0.9K are reduced, 2. 3. repetition walks, and compare root-mean-square error, directly
To Ri>Ri-1When, the depth shift factor K at this moment is optimum value.
Step 3：Using Forward Formula, resistivity transforming function transformation function T is calculatedL(k)。
Step 4：Build minimum object function δ.
The T that forward modeling is obtainedL(k) the transforming function transformation function T with surveying1(k) compare, the variance of the two is obtained.
M is T in formula1(k) hits, it is equal to ρs(r) hits, and
I is iterations, i=0,1,2 ....Model as i=0 is initial model；P is vector；N is the layer of model
Number；
If by TL(k,Pi) in Pi(initial value of ith iteration) is nearby unfolded by Tailor progression, will be omitted more than secondary
Higher order term, then
So as to obtain：
Number of parameters of the N=2n-1 for layer, Δ P in formulaiIt is j-th of parameter PjCorrection amount or correct step-length.
Step 5：According to the extreme value theory of the function of many variables, upper formula is converted, can be obtained：
Being write as matrix form is：ATA Δs P=ATΔ T, wherein：
Step 6：Above-mentioned equation is solved, just obtains meeting the minimum optimum stepsize vector Δ P of δ.Then to PiIt is corrected, asks
Obtain Pi+1Point (Pi+1=Pi+ΔP)。
Step 7：It repeats above-mentioned step to gather, iterate, Approach by inchmeal makes TL(k,Pi) really close to T1(k), so as to ask
It obtains and is more accurate to Δ P.
In inverting iteration, in order to improve iteration speed, each iteration with measured visual resistivity and is calculated regarding electricity as the following formula
The ratio between resistance rate is coefficient, adjusts the resistivity of each layer：
Specific practice is：Every time in iterative calculation, root-mean-square error R is judgediWhether ε (generally taking ε=5%) is less than,
If met the requirements, at this moment model parameter (resistivity, depth) is inverting as a result, otherwise adjusting resistivity by above formula, so
After rejudge whether have R2<ε repeats above step and gathers, until equal root error is met the requirements.
Step 8：Inversion result is obtained, draws ρ (z)~z curves.
Step 9：Ask for GDS values.GDS is asked for according to formula (1).
Step 10：Draw GDS~z curves namelyCurve.
By rightThe identification of curvilinear characteristic can more intuitively reflect the variation tendency of electrical parameter.
Step 11：Due to being limited by electric sounding data, measured line smoothing degree not enough (Fig. 5), is unfavorable for differential
Imaging.Therefore, in order to preferably hold the feature of curve consecutive variations and accurate derivation.Employ slide Smoothing Technique and
Cubic Spline Interpolation carries out data multiple structure.Then, high-density sampling is carried out to reconstructing later data, finally realized micro-
It is divided into as (such as Fig. 6,7,8,9).
Second part：High-resolution earth resistivity fast imaging method
This method is the product for being combined GDS methods with the Da Zhanuoke methods of inversion.Specific steps include：
Step 1：ρs(r)~r is seen as Da Zhanuoke curves
Step 2：It calculates
Step 3：It calculates depth z and corresponding resistivity, specific formula for calculation is as follows：
Step 4：Paint ρ (z)~z curves.
Step 5：Build the relational expression of adaptation coefficient.
For three layers of earth-electricity model, if H1And H2It is the first bed boundary of earth-electricity model and the depth of the second bed boundary respectively,
F1And F2It is that corresponding earth-electricity model carries out the first bed boundary and the second bed boundary that are obtained after Da Zhanuoke curve method invertings respectively
Depth.It analyzes through a large number of experiments, finds H1/F2And F2/F1Between, H2/F1And F1/F2Between there are certain relationship, take
X=F1/F2；For by the obtained initial results of Da Zhanuoke curve method invertings, the H of four type electrical sounding curves1/F2With
F2/F1Between functional relation be all nonlinear, it is possible to nonlinear least square fitting function lsqucurvefit come
It is fitted both the above functional relation.
Step 6：Adaptive correction coefficient exports.
For the correction coefficient of the first bed boundary,It is corrected for the second bed boundary
Coefficient, a2=x*f (x).Through with nonlinear least square fitting method to H1/F2With F2/F1And H2/F1With F1/F2Distribution function
After fitting, fitting function f (x) and interface correction coefficient a (plans of the Figure 10 for H-type curve of four kinds of earth-electricity models can be respectively obtained
Close function distribution map, as space is limited the reason of, only list H-type here).
For H-type geoelectric cross section：
First bed boundary fitting function：
First bed boundary correction coefficient：
Second bed boundary fitting function：
Second bed boundary correction coefficient：
Step 7：The ρ (z) in step 4~z curves are corrected by the correction coefficient of step 6, regenerate ρ (z)
~z curves.
Step 8：Ask for GDS values.GDS is asked for according to formula (1).
Step 9：Draw GDS~z curves namelyCurve.
By rightThe identification of curvilinear characteristic can more intuitively reflect the variation tendency of electrical parameter.
Step 10：Due to being limited by electric sounding data, measured line smoothing degree not enough (Fig. 5), is unfavorable for differential
Imaging.Therefore, in order to preferably hold the feature of curve consecutive variations and accurate derivation.Employ slide Smoothing Technique and
Cubic Spline Interpolation carries out data multiple structure.Then, high-density sampling is carried out to reconstructing later data, finally realized
GDS differential imagings (such as Figure 11,12,13,14).
Step 1~4 in method two are specifically described：
1. on first branch asymptote, ρs(r)=ρ1, a=0, therefore have ρ (z)=ρ1；
2. in ρnFor on limited tail branch asymptote, ρs(r)=ρn, a=0, so ρs(r)=ρn；
3. work as ρnDuring → ∞, the slope of tail branch asymptote is 45 ° of straight lines, at this time a=1, therefore ρ (z)=∞；
4. work as ρnWhen → 0, the slope a of tail branch asymptote<- 1 will lose meaning.Work as a during practical calculating<When -1, a should be taken
=-1.
Present disclosure is not limited to cited by embodiment, and those of ordinary skill in the art are by reading description of the invention
And to any equivalent transformation that technical solution of the present invention is taken, it is that claim of the invention is covered.
Claims (6)
1. high-resolution earth resistivity fast imaging method, it is characterised in that：
Include the following steps：
(1) inverting of high-resolution earth resistivity image reconstruction, on the basis of routinely enlightening inversion method is helped, the electricity that is obtained to inverting
Resistance rate curve is differentiated, and is drawn GDS~Z figures, is realized differential imaging；
(2) high-resolution earth resistivity fast imaging, for three layers of earth-electricity model, on the basis of the Da Zhanuoke methods of inversion, difference
Four kinds of difference Da Zhanuoke inversional curves are corrected with the fitting of coefficient, passes through nonlinear least square fitting function
1squcurvefit fits first layer and the correction coefficient of the second layer draws GDS~Z figures, realizes differential imaging.
2. high-resolution earth resistivity fast imaging method according to claim 1, it is characterised in that：
The step of inverting of high-resolution earth resistivity image reconstruction, includes：
Establish initial model；
Determine the depth of each layer；
Using Forward Formula, resistivity transforming function transformation function T is calculatedL(k)；
Build minimum object function δ；
It obtains meeting the minimum optimum stepsize vector Δ P of δ, then to PiIt is corrected, acquires Pi+1Point (Pi+1=Pi+ΔP)；
It iterates, Approach by inchmeal makes TL(k,Pi) really close to T1(k), more it is accurate to Δ P so as to acquire；
Inversion result is obtained, draws ρ (z)~z curves；
Ask for GDS values；
Draw GDS~z curves namelyCurve；
Data are carried out with multiple structure using Smoothing Technique and Cubic Spline Interpolation is slided；Then, to reconstruct later data into
Row high-density sampling, finally realizes differential imaging.
3. high-resolution earth resistivity fast imaging method according to claim 1, it is characterised in that：
The step of high-resolution earth resistivity fast imaging, includes：
ρs(r)~r is seen as Da Zhanuoke curves
It calculates
Calculate depth z and corresponding resistivity；
Paint ρ (z)~z curves；
Build the relational expression of adaptation coefficient；
Export adaptive correction coefficient；
ρ (z)~z curves are corrected by correction coefficient, regenerate ρ (z)~z curves；
Ask for GDS values；
Draw GDS~z curves namelyCurve；
Data are carried out with multiple structure using Smoothing Technique and Cubic Spline Interpolation is slided；Then, to reconstruct later data into
Row high-density sampling finally realizes GDS differential imagings.
4. high-resolution earth resistivity fast imaging method according to claim 2, it is characterised in that：
The specific steps of the inverting of high-resolution earth resistivity image reconstruction include：
Step 1：Establish initial model：
Assuming that stratum the number of plies and sounding curve points as many, every layer of resistivity is regarded for respective points on sounding curve
Resistivity, every layer of depth are equal to the electrode spacings of respective points on sounding curve multiplied by with a depth shift factor, thus
An initial model (h is obtainedi,ρi, i=1,2 ... is n)；
Step 2：Determine the depth of each layer：
Since in inversion algorithm, the depth initial value of each layer is replaced with electrode spacing multiplied by with a depth shift factor, i.e.,To accelerate inverting convergence rate in inverting, following Iteration is used to the inversion result of depth：
H(i+1)=(H(i)+ΔH(i))·K (1)
Wherein K is the depth shift factor；The selection method of K is attributed to：
1. setting initial the depth shift factor as K=1, the depth value of each layer of initial model is equal to corresponding electrode spacing value；
2. calculate the theoretical apparent resistivity value ρ of initial modelC；
3. calculate theoretical apparent resistivity and the root-mean-square error R of measured visual resistivity depth measurement data：
Wherein：I is iterations, and j is jth layer or j-th of electrode spacing, ρ0(j) it is measured visual resistivity on j-th of electrode spacing
Value, ρci(j) the theoretical apparent resistivity value to be calculated on j-th of electrode spacing of iteration j；
If 4. R0>R1, then the depth shift factor, i.e. K=0.9K are reduced, 2. 3. repetition walks, and compare root-mean-square error, until Ri>
Ri-1When, the depth shift factor K at this moment is optimum value；
Step 3：Using Forward Formula, resistivity transforming function transformation function T is calculatedL(k)；
Step 4：Build minimum object function δ；
The T that forward modeling is obtainedL(k) the transforming function transformation function T with surveying1(k) compare, the variance of the two is obtained；
M is T in formula1(k) hits, it is equal to ρs(r) hits, and
I is iterations, i=0,1,2 ...；Model as i=0 is initial model；P is vector；N is the number of plies of model；
If by TL(k,Pi) in the initial value P of ith iterationiNearby it is unfolded by Tailor progression, secondary above high-order will be omitted
, then
So as to obtain：
Number of parameters of the N=2n-1 for layer, Δ P in formulaiIt is j-th of parameter PjCorrection amount or correct step-length；
Step 5：According to the extreme value theory of the function of many variables, formula (5) is converted, can be obtained：
Being write as matrix form is：ATA Δs P=ATΔ T, wherein：
Step 6：Above-mentioned equation is solved, just obtains meeting the minimum optimum stepsize vector Δ P of δ；Then to PiIt is corrected, acquires Pi+1
Point (Pi+1=Pi+ΔP)；
Step 7：It repeats above-mentioned step to gather, iterate, Approach by inchmeal makes TL(k,Pi) really close to T1(k), so as to acquire more
It is accurate to Δ P；
In inverting iteration, in order to improve iteration speed, each iteration with measured visual resistivity and calculates apparent resistivity as the following formula
The ratio between for coefficient, adjust the resistivity of each layer：
Specific practice is：Every time in iterative calculation, root-mean-square error R is judgediWhether it is less than ε, takes ε=5%, it will if met
It asks, at this moment model parameter resistivity, depth are inverting as a result, otherwise adjusting resistivity by above formula, then rejudging is
It is no to have R2<ε repeats above step and gathers, until equal root error is met the requirements；
Step 8：Inversion result is obtained, draws ρ (z)~z curves；
Step 9：Ask for GDS values；GDS is asked for according to formula (9)：
Wherein, ziFor i-th of depth value, zi-1For (i-1)-th depth value, resistivity values of the ρ (i) for i-th of depth, ρ (i-1)
Resistivity value for (i-1)-th depth；
Step 10：Draw GDS~z curves namelyCurve；
By rightThe identification of curvilinear characteristic can more intuitively reflect the variation tendency of electrical parameter；
Step 11：Due to being limited by electric sounding data, measured line smoothing degree is inadequate, is unfavorable for differential imaging processing；
Therefore, in order to preferably hold the feature of curve consecutive variations and accurate derivation；Using slip Smoothing Technique and cubic spline letter
Number interpolation carries out data multiple structure；Then, high-density sampling is carried out to reconstructing later data, finally realizes differential imaging.
5. high-resolution earth resistivity fast imaging method according to claim 3, it is characterised in that：
The specific steps of high-resolution earth resistivity fast imaging include：
Step 1：ρs(r)~r is seen as Da Zhanuoke curves
Step 2：It calculates
Step 3：It calculates depth z and corresponding resistivity, specific formula for calculation is as follows：
Step 4：Paint ρ (z)~z curves；
Step 5：Build the relational expression of adaptation coefficient；
For three layers of earth-electricity model, if H1And H2It is the first bed boundary of earth-electricity model and the depth of the second bed boundary respectively, F1With
F2It is that corresponding earth-electricity model carries out the first bed boundary and the second bed boundary that are obtained after Da Zhanuoke curve method invertings respectively
Depth；It analyzes through a large number of experiments, finds H1/F2And F2/F1Between, H2/F1And F1/F2Between there are certain relationship, take x
=F1/F2；For by the obtained initial results of Da Zhanuoke curve method invertings, the H of four type electrical sounding curves1/F2And F2/
F1Between functional relation be all nonlinear, it is possible to intended with nonlinear least square fitting function lsqucurvefit
Close both the above functional relation；
Step 6：Adaptive correction coefficient exports；
For the correction coefficient of the first bed boundary,For the second bed boundary correction coefficient,
a2=x*f (x).Through with nonlinear least square fitting method to H1/F2With F2/F1And H2/F1With F1/F2Distribution function fitting
Afterwards, the fitting function f (x) and interface correction coefficient a of four kinds of earth-electricity models can be respectively obtained；
For H-type geoelectric cross section：
First bed boundary fitting function：
First bed boundary correction coefficient：
Second bed boundary fitting function：
f2(x)=31.0405*exp (- 43.6255*x)+4.7257*exp (- 5.5429*x)
(13)
Second bed boundary correction coefficient：
Step 7：The ρ (z) in step 4~z curves are corrected by the correction coefficient of step 6, it is bent to regenerate ρ (z)~z
Line；
Step 8：GDS values are asked for, GDS is asked for according to formula (9)；
Step 9：Draw GDS~z curves namelyCurve；
By rightThe identification of curvilinear characteristic can more intuitively reflect the variation tendency of electrical parameter；
Step 10：Due to being limited by electric sounding data, measured line smoothing degree is inadequate, is unfavorable for differential imaging processing；
Therefore, in order to preferably hold the feature of curve consecutive variations and accurate derivation；Using slip Smoothing Technique and cubic spline letter
Number interpolation carries out data multiple structure；Then, high-density sampling is carried out to reconstructing later data, finally realizes GDS differential imagings.
6. high-resolution earth resistivity fast imaging method according to claim 5, it is characterised in that：
In step 1-4：
1. on first branch asymptote, ρs(r)=ρ1, a=0, therefore have ρ (z)=ρ1；
2. in ρnFor on limited tail branch asymptote, ρs(r)=ρn, a=0, so ρs(r)=ρn；
3. work as ρnDuring → ∞, the slope of tail branch asymptote is 45 ° of straight lines, at this time a=1, therefore ρ (z)=∞；
4. work as ρnWhen → 0, the slope a of tail branch asymptote<- 1 will lose meaning；Work as a during practical calculating<When -1, a=-1 should be taken.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810038410.3A CN108169801B (en) | 2018-01-16 | 2018-01-16 | High-resolution ground resistivity rapid imaging method |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810038410.3A CN108169801B (en) | 2018-01-16 | 2018-01-16 | High-resolution ground resistivity rapid imaging method |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN108169801A true CN108169801A (en) | 2018-06-15 |
| CN108169801B CN108169801B (en) | 2020-09-15 |
Family
ID=62514798
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201810038410.3A Active CN108169801B (en) | 2018-01-16 | 2018-01-16 | High-resolution ground resistivity rapid imaging method |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN108169801B (en) |
Cited By (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110333543A (en) * | 2019-07-03 | 2019-10-15 | 山东大学 | Post non of low resistance body explanation and imaging method and system based on reflection coefficient analysis |
| CN112415602A (en) * | 2020-10-15 | 2021-02-26 | 山东大学 | Tunnel resistivity advanced detection optimization method and system based on depth resolution |
| CN112597108A (en) * | 2020-11-30 | 2021-04-02 | 核工业二0八大队 | Direct current sounding inversion data processing system and method |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US6577700B1 (en) * | 2001-06-22 | 2003-06-10 | Liang-Shih Fan | Neural network based multi-criteria optimization image reconstruction technique for imaging two- and three-phase flow systems using electrical capacitance tomography |
| CN101334484A (en) * | 2008-07-22 | 2008-12-31 | 江苏大学 | Three-dimensional high definition electric resistivity exploration and direct imaging method |
| CN102635347A (en) * | 2012-03-30 | 2012-08-15 | 中国电子科技集团公司第二十二研究所 | Method for quantitatively enabling thin interbed to be equivalent to formation with horizontal and vertical resistivities |
| CN103869173A (en) * | 2014-02-26 | 2014-06-18 | 国家电网公司 | Method for measuring earth resistivity distribution from earth surface to underground tens of kilometers |
| CN107305600A (en) * | 2016-04-21 | 2017-10-31 | 新疆维吾尔自治区煤炭科学研究所 | Least square method resistivity three-dimensional approximate inversion technology |
-
2018
- 2018-01-16 CN CN201810038410.3A patent/CN108169801B/en active Active
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US6577700B1 (en) * | 2001-06-22 | 2003-06-10 | Liang-Shih Fan | Neural network based multi-criteria optimization image reconstruction technique for imaging two- and three-phase flow systems using electrical capacitance tomography |
| CN101334484A (en) * | 2008-07-22 | 2008-12-31 | 江苏大学 | Three-dimensional high definition electric resistivity exploration and direct imaging method |
| CN102635347A (en) * | 2012-03-30 | 2012-08-15 | 中国电子科技集团公司第二十二研究所 | Method for quantitatively enabling thin interbed to be equivalent to formation with horizontal and vertical resistivities |
| CN103869173A (en) * | 2014-02-26 | 2014-06-18 | 国家电网公司 | Method for measuring earth resistivity distribution from earth surface to underground tens of kilometers |
| CN107305600A (en) * | 2016-04-21 | 2017-10-31 | 新疆维吾尔自治区煤炭科学研究所 | Least square method resistivity three-dimensional approximate inversion technology |
Cited By (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110333543A (en) * | 2019-07-03 | 2019-10-15 | 山东大学 | Post non of low resistance body explanation and imaging method and system based on reflection coefficient analysis |
| CN112415602A (en) * | 2020-10-15 | 2021-02-26 | 山东大学 | Tunnel resistivity advanced detection optimization method and system based on depth resolution |
| CN112597108A (en) * | 2020-11-30 | 2021-04-02 | 核工业二0八大队 | Direct current sounding inversion data processing system and method |
Also Published As
| Publication number | Publication date |
|---|---|
| CN108169801B (en) | 2020-09-15 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN108169801A (en) | High-resolution earth resistivity fast imaging method | |
| CN104459782B (en) | Seismic velocity modeling method and modeling unit using thin layer chromatography inversion | |
| CN106501867B (en) | A kind of transient electromagnetic inversion method based on lateral smoothness constraint | |
| CN105974479B (en) | The chromatography 2D/3D anisotropy Depth Domain velocity modeling method of GPU space lattice body | |
| CN110007357B (en) | Aviation TEM and aviation MT joint inversion method | |
| CN113552625B (en) | Multi-scale full waveform inversion method for conventional land-domain seismic data | |
| CN107884825B (en) | Uncertainty modeling method based on seismic multi-attribute | |
| CN106096081B (en) | The estimation method of reserve of fracture hole type bottom water reservoir | |
| CN102901985A (en) | Depth domain layer speed correcting method suitable for undulating surface | |
| CN102338887B (en) | Irregular-size space-variant grid tomography imaging statics correction method | |
| CN109001825A (en) | Across the hole CT of four-dimensional resistivity based on priori gradient constraint monitors imaging method | |
| BR112015000879B1 (en) | System and method for modeling migration speed | |
| CN110133716A (en) | Magnetic anomaly data three-dimensional inversion method based on built-up pattern weighting function | |
| CN104932021A (en) | Constrained tomography speed modeling method based on reverse ray tracing | |
| CN109826174A (en) | A kind of slope reinforcement deep regional range determining method | |
| WO2020087767A1 (en) | Velocity inversion method based on jointly collected station and three-dimensional seismic data | |
| CN111221035A (en) | Seismic reflection wave slope and gravity anomaly data joint inversion method | |
| CN109884700A (en) | Multi-information fusion seismic velocity modeling method | |
| CN109387868A (en) | A kind of three-dimensional chromatography imaging method based on seismic wave lineups slope information | |
| US20140095078A1 (en) | Method and system for presenting seismic information | |
| CN105223630B (en) | Omnibearing observation systematic parameter Demonstration Method based on geological model | |
| CN114236624B (en) | Method and system for estimating fracturing modification space volume based on electromagnetic method | |
| CN106846481B (en) | Geological profile generation method | |
| CN115906559A (en) | Magnetotelluric self-adaptive finite element forward modeling method based on mixed grid | |
| CN110440754B (en) | Actual measurement geological profile method based on space coordinates |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |