CN105589108B - Transient electromagnetic quick three-dimensional inversion method based on various boundary conditions - Google Patents

Transient electromagnetic quick three-dimensional inversion method based on various boundary conditions Download PDF

Info

Publication number
CN105589108B
CN105589108B CN201510926452.7A CN201510926452A CN105589108B CN 105589108 B CN105589108 B CN 105589108B CN 201510926452 A CN201510926452 A CN 201510926452A CN 105589108 B CN105589108 B CN 105589108B
Authority
CN
China
Prior art keywords
mrow
msub
infinitesimal
time constant
moment
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.)
Active
Application number
CN201510926452.7A
Other languages
Chinese (zh)
Other versions
CN105589108A (en
Inventor
武欣
饶丽婷
张晓娟
方广有
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Institute of Electronics of CAS
Original Assignee
Institute of Electronics of CAS
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Institute of Electronics of CAS filed Critical Institute of Electronics of CAS
Priority to CN201510926452.7A priority Critical patent/CN105589108B/en
Publication of CN105589108A publication Critical patent/CN105589108A/en
Application granted granted Critical
Publication of CN105589108B publication Critical patent/CN105589108B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

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)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

The invention provides a kind of transient electromagnetic quick three-dimensional inversion method based on various boundary conditions, employs the data processing based on the conversion of transient electrical magnetic moment and D integral pin-fin tube method, solves the problem that the data volume that current 3-d inversion faces is big, forward modeling is complicated；The constraints of time constant vector in inverting is constructed, inversion problem is iterated using optimization algorithm so that refutation process optimizes to the direction for being close to actual underground structure.Method proposed by the present invention can realize the quick three-dimensional inverting of subsurface anomaly target on a common computer, explain in work there is important application prospect in real-time transient electromagnetic data.

Description

Transient electromagnetic quick three-dimensional inversion method based on various boundary conditions
Technical field
The present invention relates to geophysical exploration technology underground geologic bodies detection technology field, more particularly to one kind is based on difference The transient electromagnetic quick three-dimensional inversion method of constraints.
Background technology
Transient electromagnetic method launches bi-directional pulse current using grounding electrode or earth-free loop line to underground, and underground medium exists Its inductive loop under exciting produces the secondary field changed over time, and in the intermittent phase of primary field, magnetic is measured using receiving coil Field signal, by the extraction and analysis to secondary field signal, so as to reach the purpose of Underground geologic body.Currently as main One of Non-seismic methods, be widely used in the subterranean resource field of detecting such as oil gas, mineral products.The inversion interpretation of transient electromagnetic data Work is the important step in transient electromagnetic method exploration.Due to the complexity of the positive algorithm of higher-dimension, multidimensional inversion problem is not yet appropriate Kind to solve, in practical application, the explanation work to transient electromagnetic data is concentrated mainly on one-dimensional inversion.
But higher-dimension inverting can provide more fine ground electricity structural information, at present, with the development of computing technique, state Outer to have carried out 3-d inversion work, main the methods of utilizing integral equation method and FInite Element, realizes strict three-dimensional FORWARD AND INVERSE PROBLEMS. The strict 3-d inversion method of transient electromagnetic is limited to the D integral pin-fin tube algorithm of complexity, and data volume is huge, and occupancy resource is excessive, It can not almost be run in common computer, while the arithmetic speed of strict 3-d inversion is slowly, it is necessary to even several days several hours 3-d inversion could be completed, due to these limitations, strict 3-d inversion can't really put into practical application, can not also enter in real time The explanation of row transient electromagnetic data.Therefore, the technical problem that those skilled in the art need urgently to solve is exactly：It is how real A kind of existing 3-d inversion method for putting into application, can be finally inversed by subsurface anomaly objective body efficiently and exactly.
The content of the invention
(1) technical problems to be solved
In order to solve, data volume present in existing 3-d inversion method is big, arithmetic speed is slow, is difficult to input application Technical problem, the invention provides a kind of transient electromagnetic quick three-dimensional inversion method based on various boundary conditions.
(2) technical scheme
A kind of according to an aspect of the invention, there is provided transient electromagnetic quick three-dimensional inverting based on various boundary conditions Method, including：Step A：Emitter and receiving point are laid on the ground, setting subsurface anomaly region, abnormal area are divided For the infinitesimal of cubic shaped, according to the geometric parameter computational geometry coupling factor matrix of emitter, receiving point and infinitesimal； Step B：Emitter emission current signal, after switch off current, each receiving point gathers magnetic field data, then using uniformly greatly Model, the magnetic field data of collection is converted into apparent conductivity depth map；Step C：According to transient electromagnetic first moment convert, based on regarding Electrical conductivity depth map and measurement magnetic field data obtain the reference first moment of the abnormal area at receiving point；Step D：Build in inverting The constraints of time constant vector；Step E：Time corresponding to reference first moment and its error, infinitesimal based on abnormal area The object function of constant vector and geometry coupling factor matrix structure complementary operation；Step F：It is normal based on the time corresponding to infinitesimal The constraints of number vector, the object function of complementary operation is iterated using optimization algorithm；And step G：Judgement changes Whether the object function after generation restrains, if convergence, preserves and optimizes time constant vector and be used as inverting final result, i.e. and τ= (τ1, τ2..., τK)T, otherwise, return to step F is performed, wherein, τkBe in abnormal area corresponding to the electrical conductivity of k-th of infinitesimal when Between constant, according to time constant vector inverting final result obtain position and the volume of anomalous body.
(3) beneficial effect
It can be seen from the above technical proposal that the transient electromagnetic quick three-dimensional inverting based on various boundary conditions of the present invention Method has the advantages that：
(1) data processing method based on the conversion of transient electrical magnetic moment is employed so that the data of 3-d inversion process processing Amount is greatly reduced, and accelerates arithmetic speed, and 3-d inversion is run on a common computer；
(2) simplified D integral pin-fin tube algorithm is employed, so that the computing of forward modeling part is accelerated in refutation process；
(3) optimization algorithm based on constraints is employed, makes refutation process to the side for being close to actual underground structure To optimization, detection accuracy is improved.
Brief description of the drawings
Fig. 1 is the three-dimensional geometry schematic diagram of Simulation Calculation；
Fig. 2 is transmitting wire frame and the two dimensional topology figure for receiving survey line；
Fig. 3 is the flow chart of the present embodiment quick three-dimensional inversion method；
Fig. 4 is the CDI sectional drawings of the Transient electromagnetic measure data of emulation；
Fig. 5 is measurement data first moment, background first moment, abnormal area first moment in each bar survey line upper curve figure；
Fig. 6 is the inverse time constant sectional drawing of unconfined condition；
Fig. 7 is the inverse time constant sectional drawing of CDI starting models；
Fig. 8 is the sectional drawing of electrical conductivity weighted value；
Fig. 9 is the inverse time constant sectional drawing weighted using electrical conductivity；
Figure 10 is the sectional drawing of depth weighted value；
Figure 11 is to use depth weighted inverse time constant sectional drawing；
Figure 12 is the big fixed source-loop observation device schematic diagram of rectangle.
Embodiment
For the object, technical solutions and advantages of the present invention are more clearly understood, below in conjunction with specific embodiment, and reference Accompanying drawing, the present invention is described in more detail.It should be noted that in accompanying drawing or specification description, similar or identical portion Divide and all use identical figure number.The implementation for not illustrating or describing in accompanying drawing, it is those of ordinary skill in art Known form.In addition, though the demonstration of the parameter comprising particular value can be provided herein, it is to be understood that parameter is without definite etc. In corresponding value, but can be similar to be worth accordingly in acceptable error margin or design constraint.Mentioned in embodiment Direction term, such as " on ", " under ", "front", "rear", "left", "right" etc., only it is the direction of refer to the attached drawing.Therefore, the side used Protection scope of the present invention is intended to be illustrative and not intended to limit to term.
Present invention employs a kind of transient electromagnetic quick three-dimensional inversion method based on various boundary conditions, by simplified D integral pin-fin tube algorithm, the non-linear strict inversion problem of complexity is converted into simple linear inverse problems；In order to solve directly Inverting can not provide rational result, introduce various boundary conditions, and the three-dimensional for quickly and accurately realizing subsurface anomaly body is anti- Drill.
In one exemplary embodiment of the present invention, there is provided a kind of to the Emit Maxwell with Marco modules The method that the emulation transient electromagnetic data of software implements quick three-dimensional inverting.Marco modules are by Australian Union's Scientific Industries Research institute (CSIR O) develops, and the module is based on three-dimensional integral equation, can calculate abnormal containing multiple prisms in multi-layered earth The Transient electromagnetic response of target.
Fig. 1 is the three-dimensional geometry schematic diagram of Simulation Calculation, as illustrated, being that 1mS/m is uniformly carried on the back greatly in electrical conductivity The flat board anomalous body that electrical conductivity is 1S/m is placed in scape, flat board size is 800E × 800N × 300Z (800 meters × north and south of East and West direction To 300 meters of 800 meters × depth), wherein, E, N, Z represent East and West direction, north-south, depth to the upper surface center of the flat board respectively Coordinate is (0E, 0N, -400Z).
Fig. 3 is the flow chart of the present embodiment quick three-dimensional inversion method.Fig. 3 is refer to, the present embodiment includes：
Step A, emitter and receiving point are laid on the ground, setting subsurface anomaly region, abnormal area is divided into The infinitesimal of cubic shaped, according to the geometric parameter computational geometry coupling factor matrix of emitter, receiving point and infinitesimal；
Step A is specifically included：
Sub-step A1：In detected region, emitter and receiving point are set；
Specifically, sub-step A1 is specifically included：Laying rectangle determines greatly source wire frame on the ground, and wire frame is centrally located at (0E, 0N, 0Z), the wire frame length of side are 500 meters × 500 meters, and being uniformly distributed 11 from -500N to 500N on north-south receives survey Line, survey line spacing are 100m, and the trend of every survey line is evenly distributed with 21 receiving points from -1000E to 1000E, share and receive Point N=231.
In the present embodiment, emitter is that rectangle determines greatly source transmitting wire frame.Fig. 2 is the two of transmitting wire frame and reception survey line Tie up layout.
Sub-step A2：Abnormal area is divided into infinitesimal, wherein, the abnormal area is different to there may be in detected region The target area of normal body；Specifically, sub-step A2 is specifically included：Abnormal area scope in inverting is set to arrive as -1000E 1000E (transmeridional -1000 meters to+1000 meters), -500N arrive 500N (- 500 meters of North and South direction are to+500 meters), -2000Z To -200Z (- 200 meters of depth direction are to -2000 meters), abnormal area is divided into the cube infinitesimal that the K length of side is 25m, Infinitesimal total quantity K=230400；
Sub-step A3：It is several according to the geometric parameter computational geometry coupling factor matrix of emitter, receiving point and infinitesimal What coupling factor matrix GnkFor：
In above formula,Represent that the unit direction vector of n-th of receiving point is pointed at k-th of infinitesimal center,Expression is incided The unit direction vector of the primary field of k-th of infinitesimal, Vk represent the volume of k-th of infinitesimal, B0, kK-th of infinitesimal is incided in expression Field amplitude, rnkRepresent k-th of infinitesimal to the distance of n-th of receiving point.
In the present embodiment, sub-step A3 geometry coupling factor matrix GnkSize is 231 × 230400；
In this example, the storage area that geometry coupling factor matrix takes is 390MB.The geometry coupling factor matrix and wink Change electromagnetic measurement data are unrelated, and the geometric position that it is only received with transmitting is information-related, and therefore, the matrix can be before measuring Completed with regard to calculating, as a loading data item in iterative inversion is implemented；
Step B：Emitter emission current signal, after switch off current, each receiving point gathers magnetic field data, then uses Uniform Earth model, the magnetic field data of collection is converted into apparent conductivity depth map (CDI)：
Fig. 4 is the CDI sectional drawings of the Transient electromagnetic measure data of emulation, it can be seen that high conductivity region one Determine the position that degree has reacted flat board anomalous body, but the bottom position of high conductivity region is abnormal far beyond actual plate The bottom of body.Initial time t can be obtained from CDI1With deadline tnThe apparent conductivity σ at place1And σn, and estimate background electricity Conductance σbg.The CDI can be used for calculating CDI starting models and electrical conductivity weighted value, in field survey data inversion, CDI Reference can be used as, verifies the reliability of inversion result；
Step C, converted according to transient electromagnetic first moment, obtained based on apparent conductivity depth map (CDI) and measurement magnetic field data Take the reference first moment of the abnormal area at receiving point；
Transient electromagnetic first moment transform definition is：
That is integration of the magnetic responsiveness from 0 moment to ∞.Resistive be limited in time domain of transient electromagnetic is equal to first moment, by magnetic field Measurement data is converted into resistive limitation data, equivalent first moment conversion.
Step C is specifically included：
Sub-step C1：
Measurement data first moment is calculated based on measurement data and the vertical magnetic field response of uniform big ground；
Sub-step C1 is specifically included：
Because measurement data is in limited time range, in order to obtain the first moment of measurement data, it is necessary to the polishing time Magnetic field integral part outside scope, here, the expression formula of the measurement data first moment of emulation is：
Wherein, t1And tnMeasurement data initial time and deadline, σ are represented respectively1And σnT is represented respectively1And tnPlace Apparent conductivity, head represent magnetic field from 0 to t1Integration in time, middle part represent measurement data from t1To tnIntegration, afterbody table Show magnetic field from tnTo the integration in the ∞ times.
Middle part is the fitting function for obtaining measurement data B (t) first, by fitting function in time t1And tnInterior numerical integration Middle part is produced to integrate.Because magnetic field can not measure in the time range on head and afterbody, here using in theory uniformly greatly Magnetic responsiveness integration replace, head by calculate from 0 to t1Interior magnetic responsiveness uniformly greatly integrates to obtain；Afterbody passes through meter Calculate from tnIntegrate to obtain to the magnetic responsivenesses of ∞ uniformly greatly, wherein uniformly electrical conductivity greatly is estimated from apparent conductivity depth map Meter obtains.
It is related to the specific solution of magnetic responsiveness integration in theory uniformly greatly in the measurement data first moment of above-mentioned emulation It is as follows：
G (x, y, t) is the uniform vertical magnetic field response greatly that unitary current member excites, and its expression formula is：
Wherein,σ represents electrical conductivity uniformly greatly, estimates to obtain from apparent conductivity depth map, μ is vacuum Magnetic conductivity, x and y represent position coordinates of the correspondingly received point with respect to the end of line current member,T is receiving channel Time window.
In the present embodiment, due to using the big fixed source-loop of rectangle, therefore under rectangular loop source forcing, uniform ground greatly Magnetic responsiveness Bz(t) it is：
In above formula, G (x, y, t) is that the vertical magnetic field of unitary current member responds, x1=XE-XR, x2=XW-XR, y1=YN-YR, y2=YS-YR, wherein (XE, YN) be the big fixed source-loop northeastward of rectangle apex coordinate, (XW, YS) it is the big fixed source-loop of rectangle The apex coordinate of southwestward, as shown in figure 12.
Define the integrated form of the vertical magnetic field response of unitary current member：
As t → 0, the expression formula of above formula bracket part is equal to zero, the integral part of end can further abbreviation, because This, the wherein a line frame of the big fixed source-loop of rectangle encourages the expression formula for the first moment to be formed to be
Above formula is represented from frame (x, y1) arrive (x, y2) integration, under the excitation of rectangle big fixed source-loop, uniform the one of ground greatly Rank square is equal to the first moment sum of four edges frame.
Under a line frame excitation of the big fixed source-loop of rectangle, uniformly the magnetic responsiveness on ground greatly is from the time 0 to t1Integration be：
Uniformly big magnetic responsiveness is from time tnIntegration to ∞ is：
In the present embodiment, the uniform magnetic field greatly that the big fixed source-loop each edge frame of rectangle is calculated according to expression formula (8) rings Should be from the time 0 to t1Integration, the above-mentioned integrated value of four edges frame is summed, you can try to achieve the head in measurement data first moment Portion；The uniform magnetic responsiveness greatly of the big fixed source-loop each edge frame of rectangle is calculated from time t according to expression formula (9)nArrive ∞'s Integration, the above-mentioned integrated value of four edges frame is summed, you can try to achieve the afterbody in measurement data first moment.
In order that the measurement data of emulation is close to field survey data, the measurement data first moment of above-mentioned emulation is added 5% white Gaussian noise；
Sub-step C2：The estimation background conductance rate from apparent conductivity depth map (CDI), then utilize uniform big ground perpendicular magnetic The first moment expression formula of field obtains background response first moment；
Sub-step C2 is specifically included：Initial time t is obtained from apparent conductivity depth map CDI1With deadline tnPlace regards Conductivityσ1And σn, and estimate background conductance rate σbg。
Wherein above-mentioned estimation background conductance rate σbgSpecifically include：
Background response first moment can be by uniformly the magnetic responsiveness on ground is got from 0 to ∞ time inner products greatly.
In the present embodiment, background response first moment is according to expression formula (7), by the four edges frame of the big fixed source-loop of rectangle First moment is added to obtain.
Sub-step C3：Background first moment is rejected from measurement data first moment and obtains the reference first moment of abnormal area；
Fig. 5 is measurement data first moment, background first moment, abnormal area first moment in each bar survey line upper curve figure.This is different The reference first moment in normal region is as the reference first moment d=(d in ensuing inverting1, d2..., dN)T, in the example, connect Sink number N=231.
Step D, build the constraints of time constant vector in inverting；Constraints in the step D is following three One kind in kind constraints：
Constraints a：In refutation process, the value in time constant vector be limited on the occasion of, when inverting originates, the time The initial value of constant vector is obtained by apparent conductivity depth map；
Constraints b：In refutation process, the value in time constant vector be limited on the occasion of, and according to infinitesimal regard electricity Conductance assigns weighted value to the value in time constant vector, and when inverting originates, the initial value of time constant vector is zero；
Constraints c：In refutation process, the value in time constant vector is limited on the occasion of and according to the depth of infinitesimal Weighted value is assigned to the value in time constant vector, when inverting originates, the initial value of time constant vector is zero；
Wherein,
The initial value of time constant vector in the constraints a is obtained by the following manner：
First, initial time constant vector is definedThe time constant of k-th of infinitesimal is
Wherein, σkFor the apparent conductivity value of k-th of infinitesimal, being obtained by apparent conductivity depth map CDI, μ is space permeability, L is the length of side of infinitesimal.
According to least square method, actual measurement response is with the fitting difference for originating theoretical response
Whereinβ is fitting coefficient, dnFor the reference first moment of abnormal area.To make fitting difference minimum, There are partial differential equationThe partial differential equation are solved, obtain β, initial time constant vector is τ=β τ0。
Value calculating method is weighted in the constraints b is：
Position of the electrical conductivity weighted value to infinitesimal in CDI is related, and it is distributed between 0 and 1, k-th of infinitesimal electrical conductivity Weighted value is
Wherein, σkFor the apparent conductivity value of k-th of infinitesimal, obtained by apparent conductivity depth map CDI, σmaxFor in whole CDI Maximum apparent conductivity.Electrical conductivity weighting makes big conductive regions in CDI account for leading role, weakens the geometry of small conductive regions Coupling factor influences.
The weighting value calculating method of the constraints c is：
Depth weighted value is related to the depth of infinitesimal, and the depth weighted value of k-th of infinitesimal is
W(zk, s0, d0)=tanh (s0·(zk-d0)) (13)
In above formula, zkFor the depth of k-th of infinitesimal, s0For slope factor, d0For reference depth.To different depth infinitesimal Time constant assigns weighted value, the depth weighted influence for reducing shallow-layer coupling factor.
Step E, time constant vector and geometry corresponding to reference first moment and its error, infinitesimal based on abnormal area Coupling factor matrix builds the object function of complementary operation；
Step E is specifically included：
cnThe theoretical first moment of abnormal area forward modeling computing at the n-th receiving point is represented, its computational methods is：
In theory, the total first moment response of the earth is made up of background response first moment and anomalous body response first moment superposition, table It is as follows up to formula：
In above formula,The background response first moment at n-th of receiving point is represented, part of being summed on the right of equation represents abnormal Body responds first moment, GnkGeometry coupling factor matrix for k-th of infinitesimal relative to n-th of receiving point, τkFor k-th infinitesimal Time constant.Total first moment of the earth should be equal to measurement data first moment in theory, herein, need to be from measurement data first moment It is finally inversed by unknown τk。
For simplified operation, background response first moment is rejected from measurement data first moment in refutation process, exceptions area The theoretical first moment of domain forward modeling computing takes：
Then the object function expression formula of complementary operation is
Wherein, the time constant vector of infinitesimal is τ=(τ1, τ2..., τK)T, τkIt is abnormal for the time constant of k-th of infinitesimal The error of the reference first moment in region is q=(q1, q2..., qN)T, qnFor the error of the reference first moment of n-th of abnormal area, cnRepresent the theoretical first moment of abnormal area forward modeling computing at the n-th receiving point, GnkFor geometry coupling factor matrix, dnFor exceptions area The reference first moment in domain.
Step F, based on the constraints of time constant vector corresponding to infinitesimal, using optimization algorithm to complementary operation Object function is iterated；
Step F is specifically included：
Sub-step F1：
According to constraints a, then the time constant vector before i+1 time iteration is τi, the time after i+1 time iteration Constant vector is τi+1=τi+δτi；
According to constraints b, then the time constant vector before i+1 time iteration is τi, the time after i+1 time iteration Constant vector is τi+1=τi+W(σk)δτi；
According to constraints c, then the time constant vector before i+1 time iteration is τi, the time after i+1 time iteration Constant vector is τi+1=τi+W(zk, s0, d0)δτi；
Sub-step F2：
The object function is iterated using optimization algorithm.
Optimization algorithm is one kind in following methods：Steepest descent method, conjugate gradient method, Newton method, quasi-Newton method, letter Rely domain method, singular value decomposition method, simulated annealing, genetic algorithm.
In the present embodiment, alternative manner is illustrated by taking steepest descent method as an example,
δτiCalculated by following formula：
Step-length αiExpression formula be
In above formula, intermediate variableWithExpression formula be
Wherein,It is the gradient for being fitted poor object function, is expressed as
Step G, judges whether the object function after iteration restrains, if convergence, preserves and optimizes time constant vector work For inverting final result, i.e. τ=(τ1, τ2..., τK)T, otherwise, return to step F is performed, wherein, τkIt is k-th in abnormal area The corresponding time constant of the electrical conductivity of infinitesimal, the position of anomalous body is obtained according to the inverting final result of time constant vector.
Wherein, whether the object function, which restrains, refers to whether object function is less than 1, less than 1 convergence, does not otherwise receive Hold back.
TimeconstantτkRelevant with the electrical conductivity of corresponding k-th of infinitesimal, electrical conductivity is bigger, and time constant is bigger.
Fig. 6 is the inverse time constant sectional drawing of unconfined condition, and the inversion result passes through 73 iteration, takes 60s.By There is negative value in the unconfined condition in refutation process, inverse time constant, meanwhile, because the geometry coupling factor of shallow-layer is bigger, Abnormal area ground proximity in inversion result.It can be seen that nonrestraint inversion PSD can not reflect the truly electric structure of flat board anomalous body, Therefore, in refutation process, setting constraints is very important.
Fig. 7 is the inverse time constant sectional drawing for the CDI starting models that the constraints a in step D is related to.
The initial fitting difference of CDI starting models is 1.58, iteration 8 times, and fitting difference is reached after time 0.5s less than the receipts for 1 Hold back standard.From figure 7 it can be seen that time constant inversion result meeting convergence CDI forms, the time constant maximum and flat board of inverting Actual time constant value (the τ of anomalous bodyslab=4.5msec) it coincide very much.
Fig. 8 is the sectional drawing for the electrical conductivity weighted value that the constraints b in step D is related to.
Electrical conductivity weighting makes big conductive regions in CDI account for leading role, weaken the geometry couplings of small conductive regions because Son influences.
Fig. 9 is the inverse time constant sectional drawing weighted using electrical conductivity.Initial fitting difference is 67, after inverting iteration 50 times Reach convergence, take 43s, as illustrated, different from the inversion result of CDI starting models, electrical conductivity weights high in inversion result Time constant region is confined in the volume range of actual plate, accurately reflects position and the volume of anomalous body.
Figure 10 is the sectional drawing for the depth weighted value that the constraints c in step D is related to.
In this example, s0=0.004, d0=350, the depth weighted value that the setup parameter calculates is as shown in Figure 10.
Figure 11 is to use depth weighted inverse time constant sectional drawing.Initial fitting difference is 67, reaches receipts after iteration 5 times Hold back, take 8s, as illustrated, high time constant region preferably coincide with the position of flat board anomalous body and volume.
So far, the present embodiment is described in detail combined accompanying drawing.According to above description, those skilled in the art There should be clear understanding to transient electromagnetic quick three-dimensional inversion method of the present invention based on various boundary conditions.
In addition, the above-mentioned definition to each element and method is not limited in various concrete structures, the shape mentioned in embodiment Shape or mode, those of ordinary skill in the art simply can be changed or replaced to it.
In summary, the invention provides a kind of transient electromagnetic quick three-dimensional inversion method based on various boundary conditions, The process employs a kind of processing method of data compression and the D integral pin-fin tube algorithm of simplification, solve current 3-d inversion and face Data volume is big, problem that forward modeling is complicated；Less qualitative due to 3-d inversion problem, direct inversion can not typically provide reasonably As a result,, can be more accurate method proposes different constraints in order to promote inversion result closing to reality underground structure It is finally inversed by subsurface anomaly body.Method proposed by the present invention can realize that the quick three-dimensional of subsurface anomaly target is anti-on a common computer Drill, explain in work that there is important application prospect in real-time transient electromagnetic data.
Particular embodiments described above, the purpose of the present invention, technical scheme and beneficial effect are carried out further in detail Describe in detail it is bright, should be understood that the foregoing is only the present invention specific embodiment, be not intended to limit the invention, it is all Within the spirit and principles in the present invention, any modification, equivalent substitution and improvements done etc., it should be included in the guarantor of the present invention Within the scope of shield.

Claims (9)

  1. A kind of 1. transient electromagnetic quick three-dimensional inversion method based on various boundary conditions, it is characterised in that including：
    Step A：Lay emitter and receiving point on the ground, setting subsurface anomaly region, abnormal area is divided into cube The infinitesimal of shape, according to the geometric parameter computational geometry coupling factor matrix of emitter, receiving point and infinitesimal；It is described several What coupling factor matrix GnkFor：
    <mrow> <msub> <mi>G</mi> <mrow> <mi>n</mi> <mi>k</mi> </mrow> </msub> <mo>=</mo> <msub> <mi>B</mi> <mrow> <mn>0</mn> <mo>,</mo> <mi>k</mi> </mrow> </msub> <msub> <mi>V</mi> <mi>k</mi> </msub> <mfrac> <mrow> <mn>3</mn> <mrow> <mo>(</mo> <msub> <mover> <mi>b</mi> <mo>^</mo> </mover> <mi>k</mi> </msub> <mo>&amp;times;</mo> <msub> <mover> <mi>r</mi> <mo>^</mo> </mover> <mrow> <mi>n</mi> <mi>k</mi> </mrow> </msub> <mo>)</mo> </mrow> <msub> <mover> <mi>r</mi> <mo>^</mo> </mover> <mrow> <mi>n</mi> <mi>k</mi> </mrow> </msub> <mo>-</mo> <msub> <mover> <mi>b</mi> <mo>^</mo> </mover> <mi>k</mi> </msub> </mrow> <mrow> <mn>4</mn> <msup> <msub> <mi>&amp;pi;r</mi> <mrow> <mi>n</mi> <mi>k</mi> </mrow> </msub> <mn>3</mn> </msup> </mrow> </mfrac> <mfrac> <msup> <mi>&amp;pi;</mi> <mn>2</mn> </msup> <mn>10</mn> </mfrac> </mrow>
    In above formula,Represent that the unit direction vector of n-th of receiving point is pointed at k-th of infinitesimal center,Expression is incided k-th The unit direction vector of the primary field of infinitesimal, VkRepresent the volume of k-th of infinitesimal, B0,kK-th of infinitesimal primary field is incided in expression Amplitude, rnkRepresent k-th of infinitesimal to the distance of n-th of receiving point；
    Step B：Emitter emission current signal, after switch off current, each receiving point gathers magnetic field data, then using uniform Earth model, the magnetic field data of collection is converted into apparent conductivity depth map；
    Step C：Converted according to transient electromagnetic first moment, obtained based on apparent conductivity depth map and measurement magnetic field data at receiving point Abnormal area reference first moment；
    Step D：Build the constraints of time constant vector in inverting；
    Step E：Time constant vector and geometry coupling corresponding to reference first moment and its error, infinitesimal based on abnormal area Factor matrix builds the object function of complementary operation；
    Step F：Based on the constraints of time constant vector corresponding to infinitesimal, the target using optimization algorithm to complementary operation Function is iterated；And
    Step G：Judge whether the object function after iteration restrains, if convergence, preserve and optimize time constant vector as anti- Drill final result, i.e. τ=(τ1,τ2,…,τK)T, otherwise, return to step F is performed, wherein, τkIt is k-th of infinitesimal in abnormal area Electrical conductivity corresponding to time constant, according to time constant vector inverting final result obtain position and the volume of anomalous body.
  2. 2. transient electromagnetic quick three-dimensional inversion method according to claim 1, it is characterised in that the pact in the step D Beam condition is one kind in following three kinds of constraints：
    Constraints a：In refutation process, the value in time constant vector be limited on the occasion of, when inverting originates, time constant The initial value of vector is obtained by apparent conductivity depth map；
    Constraints b：In refutation process, the value in time constant vector is limited on the occasion of and according to the apparent conductivity of infinitesimal Weighted value is assigned to the value in time constant vector, when inverting originates, the initial value of time constant vector is zero；
    Constraints c：In refutation process, the value in time constant vector be limited on the occasion of, and according to the depth of infinitesimal to when Between value in constant vector assign weighted value, when inverting originates, the initial value of time constant vector is zero.
  3. 3. transient electromagnetic quick three-dimensional inversion method according to claim 2, it is characterised in that in the constraints a Time constant vector initial value obtained by the following manner
    Define initial time constant vectorThe time constant of k-th of infinitesimal is：
    <mrow> <msubsup> <mi>&amp;tau;</mi> <mi>k</mi> <mn>0</mn> </msubsup> <mo>=</mo> <mfrac> <mrow> <msub> <mi>&amp;sigma;</mi> <mi>k</mi> </msub> <msup> <mi>&amp;mu;L</mi> <mn>2</mn> </msup> </mrow> <mn>25.6</mn> </mfrac> </mrow>
    Wherein, σkFor the apparent conductivity value of k-th of infinitesimal, obtained by apparent conductivity depth map, μ is space permeability, and L is infinitesimal The length of side；
    According to least square method, actual measurement response is with the fitting difference for originating theoretical response：
    <mrow> <mi>S</mi> <mo>=</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>n</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>N</mi> </munderover> <msup> <mrow> <mo>(</mo> <msub> <mi>d</mi> <mi>n</mi> </msub> <mo>-</mo> <mi>&amp;beta;</mi> <mo>&amp;CenterDot;</mo> <msubsup> <mi>c</mi> <mi>n</mi> <mn>0</mn> </msubsup> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mrow>
    Whereinβ is fitting coefficient, dnFor the reference first moment of abnormal area, GnkFor geometry coupling factor square Battle array, to make fitting difference minimum, there are partial differential equationThe partial differential equation are solved, obtain β, initial time constant vector For τ=β τ0。
  4. 4. transient electromagnetic quick three-dimensional inversion method according to claim 2, it is characterised in that in the constraints b Weighting value calculating method is：
    Position of the electrical conductivity weighted value to infinitesimal in apparent conductivity depth map is related, and it is distributed between 0 and 1, k-th of infinitesimal Electrical conductivity weighted value is：
    <mrow> <mi>W</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;sigma;</mi> <mi>k</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <msub> <mi>&amp;sigma;</mi> <mi>k</mi> </msub> <msub> <mi>&amp;sigma;</mi> <mrow> <mi>m</mi> <mi>a</mi> <mi>x</mi> </mrow> </msub> </mfrac> </mrow>
    Wherein, σkFor the apparent conductivity value of k-th of infinitesimal, obtained by apparent conductivity depth map, σmaxFor whole apparent conductivity depth Maximum apparent conductivity in figure.
  5. 5. transient electromagnetic quick three-dimensional inversion method according to claim 2, it is characterised in that the constraints c's Weighting value calculating method is：
    Depth weighted value is related to the depth of infinitesimal, and the depth weighted value of k-th of infinitesimal is：
    W(zk,s0,d0)=tanh (s0·(zk-d0))
    In above formula, zkFor the depth of k-th of infinitesimal, s0For slope factor, d0For reference depth.
  6. 6. transient electromagnetic quick three-dimensional inversion method according to claim 2, it is characterised in that the step F is specifically wrapped Include：
    Sub-step F1：
    According to constraints a, then the time constant vector before i+1 time iteration is τi, the time constant after i+1 time iteration Vector is τi+1=τi+δτi；
    According to constraints b, then the time constant vector before i+1 time iteration is τi, the time constant after i+1 time iteration Vector is τi+1=τi+W(σk)δτi；
    According to constraints c, then the time constant vector before i+1 time iteration is τi, the time constant after i+1 time iteration Vector is τi+1=τi+W(zk,s0,d0)δτi；And
    Sub-step F2：
    The object function is iterated using optimization algorithm.
  7. 7. the transient electromagnetic quick three-dimensional inversion method according to any one of claim 1 to 6 claim, its feature exist In the step A is specifically included：
    Sub-step A1：In detected region, emitter and receiving point are set；
    Sub-step A2：Abnormal area is divided into infinitesimal, wherein, the abnormal area is that there may be anomalous body in detected region Target area；And
    Sub-step A3：According to the geometric parameter computational geometry coupling factor matrix of emitter, receiving point and infinitesimal.
  8. 8. transient electromagnetic quick three-dimensional inversion method according to claim 7, it is characterised in that the step E is specifically wrapped Include：
    The object function expression formula of complementary operation is：
    <mrow> <mi>S</mi> <mrow> <mo>(</mo> <mi>&amp;tau;</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mn>1</mn> <mi>N</mi> </mfrac> <munderover> <mi>&amp;Sigma;</mi> <mrow> <mi>n</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>N</mi> </munderover> <msup> <mrow> <mo>(</mo> <mfrac> <mrow> <msub> <mi>d</mi> <mi>n</mi> </msub> <mo>-</mo> <msub> <mi>c</mi> <mi>n</mi> </msub> </mrow> <msub> <mi>q</mi> <mi>n</mi> </msub> </mfrac> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mrow>
    Wherein, the time constant vector of infinitesimal is τ=(τ1,τ2,…,τK)T, τkFor the time constant of k-th of infinitesimal, abnormal area The error of reference first moment be q=(q1,q2,…,qN)T, qnFor the error of the reference first moment of n-th of abnormal area, cnTable Show the theoretical first moment of abnormal area forward modeling computing at the n-th receiving point,GnkFor geometry coupling factor matrix, dn For the reference first moment of abnormal area.
  9. 9. the transient electromagnetic quick three-dimensional inversion method according to any one of claim 1 to 6 claim, its feature exist In the step C is specifically included：
    Sub-step C1：Measurement data first moment is calculated based on measurement data and the vertical magnetic field response of uniform big ground；
    Sub-step C2：Background conductance rate is estimated from apparent conductivity depth map, then utilizes the single order of uniform big ground vertical magnetic field Square expression formula obtains background response first moment；And
    Sub-step C3：Background first moment is rejected from measurement data first moment and obtains the reference first moment of abnormal area.
CN201510926452.7A 2015-12-14 2015-12-14 Transient electromagnetic quick three-dimensional inversion method based on various boundary conditions Active CN105589108B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201510926452.7A CN105589108B (en) 2015-12-14 2015-12-14 Transient electromagnetic quick three-dimensional inversion method based on various boundary conditions

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201510926452.7A CN105589108B (en) 2015-12-14 2015-12-14 Transient electromagnetic quick three-dimensional inversion method based on various boundary conditions

Publications (2)

Publication Number Publication Date
CN105589108A CN105589108A (en) 2016-05-18
CN105589108B true CN105589108B (en) 2017-11-21

Family

ID=55928832

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201510926452.7A Active CN105589108B (en) 2015-12-14 2015-12-14 Transient electromagnetic quick three-dimensional inversion method based on various boundary conditions

Country Status (1)

Country Link
CN (1) CN105589108B (en)

Families Citing this family (13)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106646625A (en) * 2016-09-27 2017-05-10 中国科学院电子学研究所 Transient electromagnetic inversion method for sharp boundary model
CN108802834B (en) * 2018-02-13 2020-12-22 中国科学院电子学研究所 Underground target identification method based on joint inversion
CN110222429A (en) * 2019-06-10 2019-09-10 清华大学 The optimization method of more line current Reconstructions
CN110865414B (en) * 2019-11-01 2021-02-02 吉林大学 Transient electromagnetic noise suppression method for urban underground space detection
CN110989002B (en) * 2019-11-07 2021-01-05 吉林大学 Shallow low-resistance abnormal body data interpretation method for earth-space time domain electromagnetic system
CN111796335B (en) * 2020-08-28 2023-03-14 核工业航测遥感中心 Aviation transient electromagnetic time constant extraction method
CN112965120B (en) * 2021-02-01 2024-04-12 西双版纳景海高速公路建设投资有限公司 Transient electromagnetic signal processing method and device for geological exploration and storage medium
CN112949134B (en) * 2021-03-09 2022-06-14 吉林大学 Earth-well transient electromagnetic inversion method based on non-structural finite element method
CN113176617A (en) * 2021-03-15 2021-07-27 中煤科工集团西安研究院有限公司 Sedimentary stratum transient electromagnetic multi-parameter constraint inversion imaging method
CN113552637B (en) * 2021-07-30 2023-11-14 中国自然资源航空物探遥感中心 Collaborative three-dimensional inversion method for magnetic anomaly data in aviation-ground-well
CN113391362B (en) * 2021-08-13 2021-10-29 成都理工大学 Magnetotelluric profile three-dimensional structured inversion method based on corridor data constraint
CN114114438B (en) * 2021-09-27 2023-07-18 中国地质科学院地球物理地球化学勘查研究所 Quasi-three-dimensional inversion method for ground-to-air transient electromagnetic data of loop source
CN114265124B (en) * 2021-12-30 2022-08-02 成都理工大学 Unfavorable geologic body positioning method based on time domain transient electromagnetic probability inversion

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103760614A (en) * 2014-02-24 2014-04-30 中国科学院电子学研究所 Transient electromagnetic forward modeling method applicable to irregular transmitted waveforms
CN103777248A (en) * 2014-02-08 2014-05-07 中国科学院电子学研究所 TEM one-dimensional forward modeling method applicable to irregular transmitting loop
CN104360404A (en) * 2014-11-27 2015-02-18 中国科学院电子学研究所 Magnetotelluric regularization inversion method based on different constraint conditions

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8095345B2 (en) * 2009-01-20 2012-01-10 Chevron U.S.A. Inc Stochastic inversion of geophysical data for estimating earth model parameters
CN102419454A (en) * 2011-06-30 2012-04-18 中国科学院地质与地球物理研究所 Method for carrying out transient electromagnetic forecasting on long-distance water-containing target body in front of tunnel face
CN102901989A (en) * 2011-07-29 2013-01-30 中国国土资源航空物探遥感中心 Gravity field or magnetic field data based geologic body three-dimensional visualized modeling and interpretation method

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103777248A (en) * 2014-02-08 2014-05-07 中国科学院电子学研究所 TEM one-dimensional forward modeling method applicable to irregular transmitting loop
CN103760614A (en) * 2014-02-24 2014-04-30 中国科学院电子学研究所 Transient electromagnetic forward modeling method applicable to irregular transmitted waveforms
CN104360404A (en) * 2014-11-27 2015-02-18 中国科学院电子学研究所 Magnetotelluric regularization inversion method based on different constraint conditions

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
时间域电磁勘探数据的模拟退火法反演研究;张建国，等;《电子与信息学报》;20150131;第37卷(第1期);220-225 *

Also Published As

Publication number Publication date
CN105589108A (en) 2016-05-18

Similar Documents

Publication Publication Date Title
CN105589108B (en) Transient electromagnetic quick three-dimensional inversion method based on various boundary conditions
Yang et al. Three-dimensional inversion of airborne time-domain electromagnetic data with applications to a porphyry deposit
Grayver et al. 3D inversion and resolution analysis of land-based CSEM data from the Ketzin CO 2 storage formation
CN101609169B (en) Method for improving electromagnetic wave resistivity measurement precision and expanding measurement range thereof
CN112949134B (en) Earth-well transient electromagnetic inversion method based on non-structural finite element method
CN105319589B (en) A kind of fully automatic stereo chromatography conversion method using local lineups slope
CN104360401B (en) A kind of transient electromagnetic B field descends objective body geological information method definitely
CN106291725B (en) A kind of method of fast inversion underground geologic bodies spatial position
CN105386756B (en) A method of brittle formation porosity is calculated using dependent variable
CN106199742A (en) A kind of Frequency-domain AEM 2.5 ties up band landform inversion method
WO2009146041A1 (en) Constructing a reduced order model of an electromagnetic response in a subterranean structure
CN103064124B (en) A kind of ratio approach correcting electromagnetic survey Considering Terrain Effect
CN104854480A (en) Apparatus and methods to find a position in an underground formation
CN103777248A (en) TEM one-dimensional forward modeling method applicable to irregular transmitting loop
CN104408021A (en) Electric dipole source three-dimensional time domain finite difference direct interpretation imaging method
CN106483559A (en) A kind of construction method of subsurface velocity model
CN115238550A (en) Self-adaptive unstructured grid landslide rainfall geoelectric field numerical simulation calculation method
CN110458129A (en) Nonmetallic mine recognition methods based on depth convolutional neural networks
CN105573963B (en) A kind of horizontal uneven texture reconstructing method in ionosphere
CN104267443B (en) Magnetotelluric field static displacement correction method based on inversion model
CN104316961A (en) Method for obtaining geological parameters of weathered layer
Zhang et al. 3D inversion of large-scale frequency-domain airborne electromagnetic data using unstructured local mesh
CN104020508A (en) Direct ray tracking algorithm for geological radar chromatography detecting
Geng et al. 3D joint inversion of gravity-gradient and borehole gravity data
CN105550442B (en) Data processing and D integral pin-fin tube method based on the transformation of transient electrical magnetic moment

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant