CN105589108A - Rapid three-dimensional inversion method for transient electromagnetism based on different constraint conditions - Google Patents

Abstract

The invention provides a rapid three-dimensional inversion method for transient electromagnetism based on different constraint conditions. According to the technical scheme of the invention, the data processing and three-dimensional forward calculation method based on the transient electromagnetic moment inversion is adopted. In this way, the problem that the current three-dimensional inversion is large in data volume and complex in forward calculation can be solved. According to the method, firstly, a constraint condition for a time constant vector during the inversion process is constructed. After that, the iteration on the inversion process is conducted based on the optimization algorithm, so that the inversion process is optimized to be close to an actual underground structure. According to the technical scheme of the invention, based on the above method, the rapid three-dimensional inversion process of an underground abnormal target can be realized by an ordinary computer. Therefore, the method has an important application prospect during the interpretation process of 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 Techniques field, relate in particular to onePlant the transient electromagnetic quick three-dimensional inversion method based on various boundary conditions.

Background technology

Transient electromagnetic method utilizes earth electrode or earth-free loop line to underground transmitting bi-directional pulse current,The inductive loop of underground medium under it excites produces time dependent secondary field, between primary fieldHave a rest the phase, use receiving coil to measure field signal, by the extraction to secondary field signal and analysis, fromAnd reach the object of Underground geologic body. As one of main Non-seismic methods, extensively should at presentFor subterranean resource field of detecting such as oil gas, mineral products. The inversion interpretation work of transient electromagnetic data is winksImportant step in power transformation magnetic prospecting. Due to the complexity of the positive algorithm of higher-dimension, multidimensional inversion problemNot yet properly settle, in practical application, the interpretation work of transient electromagnetic data is mainly concentrated on to one dimensionInverting.

But higher-dimension inverting can provide more meticulous ground electricity structural information, at present, along with calculating skillThe development of art, has carried out 3-d inversion work abroad, mainly utilizes integral equation method and FInite ElementRealize strict three-dimensional FORWARD AND INVERSE PROBLEMS etc. method. The strict 3-d inversion method of transient electromagnetic is limited to complicatedThree-dimensional positive algorithm, data volume is huge, takies resource too much, in common computer, almost cannot transportOK, the arithmetic speed of strict 3-d inversion is slow simultaneously, needs even within several days, just can complete for several hours3-d inversion, due to these restriction, strict 3-d inversion can't really drop into practical application, also withoutMethod is carried out the explanation of transient electromagnetic data in real time. Therefore, those skilled in the art need urgent solutionA technical problem is exactly: how to realize a kind of 3-d inversion method that drops into application, and can be efficientAnd be finally inversed by exactly subsurface anomaly objective body.

Summary of the invention

(1) technical problem that will solve

In order to solve, the data volume existing in existing 3-d inversion method is large, arithmetic speed is slow, be difficult toThe technical problem that drops into application, the invention provides a kind of transient electromagnetic based on various boundary conditions fastSpeed 3-d inversion method.

(2) technical scheme

According to an aspect of the present invention, provide a kind of transient electromagnetic based on various boundary conditions fastSpeed 3-d inversion method, comprising: steps A: lay on the ground emitter and acceptance point, setSubsurface anomaly region, is divided into abnormal area the infinitesimal of cube shaped, according to emitter, connectThe geometric parameter computational geometry coupling factor matrix of sink and infinitesimal; Step B: emitter transmittingCurrent signal, electric current closes has no progeny, and each acceptance point gathers magnetic field data, then adopts evenly large topotypeType, is converted to apparent conductivity depth map by the magnetic field data of collection; Step C: according to transient electromagnetic oneRank square converts, and obtains the abnormal area at acceptance point place based on apparent conductivity depth map and measurement magnetic field dataReference first moment; Step D: the constraints that builds time constant vector in inverting; Step e:The time constant vector that reference first moment based on abnormal area and error thereof, infinitesimal are corresponding and how muchCoupling factor matrix builds the object function of complementary operation; Step F: the time based on infinitesimal is corresponding is normalThe constraints of number vector, adopts optimization algorithm to carry out iteration to the object function of complementary operation; WithAnd step G: judge whether the object function after iteration restrains, if the optimization time is preserved in convergenceConstant vector is as inverting final result, i.e. τ=(τ₁，τ₂，…，τ_K)^T, otherwise, return to step F and carry out,Wherein, τ_kBe the time constant corresponding to electrical conductivity of k infinitesimal in abnormal area, normal according to the timeThe inverting final result of number vector obtains position and the volume of anomalous body.

(3) beneficial effect

Can find out from technique scheme, the transient electromagnetic based on various boundary conditions of the present invention is fastSpeed 3-d inversion method has following beneficial effect:

(1) adopt the data processing method based on the conversion of transient electrical magnetic moment, made 3-d inversion processThe data volume of processing significantly reduces, and has accelerated arithmetic speed, and makes the 3-d inversion can be at common meterOn calculation machine, move;

(2) adopted the positive algorithm of three-dimensional of simplifying, thereby made just drilling in refutation process the computing of partAccelerate;

(3) adopted the optimization algorithm based on constraints, made refutation process to being close to practicallyThe direction optimization of lower structure, has improved detection accuracy.

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 that receives survey line;

Fig. 3 is the flow chart of the present embodiment quick three-dimensional inversion method;

Fig. 4 is the CDI sectional drawing of the transient electromagnetic measurement data of emulation;

Fig. 5 is that measurement data first moment, background first moment, abnormal area first moment are on each survey lineCurve map;

Fig. 6 is the inverting time constant sectional drawing of unconfined condition;

Fig. 7 is the inverting time constant sectional drawing of the initial model of CDI;

Fig. 8 is the sectional drawing of electrical conductivity weighted value;

Fig. 9 is the inverting time constant sectional drawing that adopts electrical conductivity weighting;

Figure 10 is the sectional drawing of depth weighted value;

Figure 11 adopts depth weighted inverting time constant sectional drawing;

Figure 12 is that rectangle is determined greatly source loop line observation device schematic diagram.

Detailed description of the invention

For making the object, technical solutions and advantages of the present invention clearer, below in conjunction with concrete realityExecute example, and with reference to accompanying drawing, the present invention is described in more detail. It should be noted that, at accompanying drawing orDuring description is described, similar or identical part is all used identical figure number. In accompanying drawing, do not illustrate or retouchThe implementation of stating is form known to a person of ordinary skill in the art in affiliated technical field. In addition, althoughCan provide the demonstration of the parameter that comprises particular value herein, but should be appreciated that, parameter is without definitely equaling correspondingValue, but can in acceptable error margin or design constraint, be similar to corresponding value. EmbodimentIn the direction term mentioned, for example " on ", D score, 'fornt', 'back', " left side ", " right side " etc., onlyIt is the direction with reference to accompanying drawing. Therefore, the direction term of use is to be not used for limiting this for explanationBright protection domain.

The present invention has adopted a kind of transient electromagnetic quick three-dimensional inversion method based on various boundary conditions,By the positive algorithm of three-dimensional of simplifying, convert the non-linear strict inversion problem of complexity to simple lineProperty inversion problem; Can not provide rational result in order to solve direct inversion, introduce various boundary conditions,Realize quickly and accurately the 3-d inversion of subsurface anomaly body.

In one exemplary embodiment of the present invention, provide a kind of to Marco moduleThe emulation transient electromagnetic data of EmitMaxwell software are implemented the method for quick three-dimensional inverting. MarcoModule is by Science Industry Research Inst of Australian Union (O of CSIR) exploitation, and this module is long-pending based on three-dimensionalDivide equation, can calculate the Transient electromagnetic response that contains multiple prism abnormal objects in stratiform the earth.

Fig. 1 is the three-dimensional geometry schematic diagram of Simulation Calculation, as shown in the figure, is 1mS/m in electrical conductivityEvenly in background, place the dull and stereotyped anomalous body that electrical conductivity is 1S/m greatly, dull and stereotyped size is800E × 800N × 300Z (300 meters of the 800 meters × degree of depth in 800 meters × north-south of East and West direction), wherein, E,N, Z represent respectively East and West direction, north-south, the degree of depth to, this dull and stereotyped upper surface centre coordinate is(0E，0N，-400Z)。

Fig. 3 is the flow chart of the present embodiment quick three-dimensional inversion method. Please refer to Fig. 3, the present embodimentComprise:

Steps A, lays emitter and acceptance point on the ground, sets subsurface anomaly region, by differentNormal region is divided into the infinitesimal of cube shaped, according to how much of emitter, acceptance point and infinitesimalHow much coupling factor matrixes of calculation of parameter;

Steps A specifically comprises:

Sub-step A1: be detected region division emitter and acceptance point;

Particularly, this sub-step A1 specifically comprises: lay on the ground rectangle and determine greatly source wire frame,Wire frame is centered close to (0E, 0N, 0Z), and the wire frame length of side is 500 meters × 500 meters, on north-south from-500NTo 500N be uniformly distributed 11 receive surveys line, interval of survey line is 100m, the trend of every survey line from-1000E, to 1000E, is evenly distributed with 21 acceptance points, total acceptance point N=231.

In the present embodiment, emitter is that rectangle is determined greatly source transmitting wire frame. Fig. 2 for transmitting wire frame andReceive the two dimensional topology figure of survey line.

Sub-step A2: abnormal area is divided into infinitesimal, and wherein, this abnormal area is for being detected districtIn territory, may there is the target area of anomalous body; Particularly, this sub-step A2 specifically comprises: establishDetermine in inverting abnormal area scope for-1000E is to 1000E (transmeridional-1000 meter extremely+1000 meters),-500N is to 500N (North and South direction-500 meters to+500 meters), and-2000Z is to-200Z (degree of depth sideTo-200 meters to-2000 meters), abnormal area is divided into the cube that K the length of side is 25m micro-Unit, infinitesimal total quantity K=230400;

Sub-step A3: according to the geometric parameter computational geometry coupling of emitter, acceptance point and infinitesimalClose factor matrix, how much coupling factor matrix G_nkFor:

$G_{nk}=B_{0,k}{V}_k\frac{3({\hat{b}}_k\×{\hat{r}}_{nk}){\hat{r}}_{nk}-{\hat{b}}_k}{4{{\mathrm{\πr}}_{nk}}^3}\frac{{\mathrm{\π}}^2}{10}---\left(1\right)$

In above formula,Represent that k infinitesimal center point to the unit direction vector of n acceptance point,The unit direction vector that represents the primary field that incides k infinitesimal, Vk represents the body of k infinitesimalLong-pending, B_0，kRepresent to incide k infinitesimal primary field amplitude, r_nkRepresent that k infinitesimal is to nThe distance of acceptance point.

In the present embodiment, how much coupling factor matrix G of sub-step A3_nkSize is 231 × 230400;

In this example, the storage area that geometry coupling factor matrix takies is 390MB. This geometry couplingFactor matrix and transient electromagnetic measurement data are irrelevant, and it is only information-related with the geometric position that transmitting receives,Therefore, this matrix can just calculate before measuring, and in enforcement iterative inversion, added as oneCarry data item;

Step B: emitter emission current signal, electric current closes has no progeny, and each acceptance point gathers magnetic fieldData, then adopt even Earth model, and the magnetic field data of collection is converted to apparent conductivity depth map(CDI)：

Fig. 4 is the CDI sectional drawing of the transient electromagnetic measurement data of emulation, as can be seen from the figure, and heightThe position of dull and stereotyped anomalous body has to a certain degree been reacted in electrical conductivity region, but the bottom in high conductivity regionPosition is far beyond the bottom of actual plate anomalous body. From CDI, can obtain initial time t₁WithDeadline t_nThe apparent conductivity σ at place₁And σ_n, and estimate background conductance rate σ_bg. This CDI can useIn calculating the initial model of CDI and electrical conductivity weighted value, in field survey data inversion, CDI also canWith as a reference, verify the reliability of inversion result;

Step C, according to transient electromagnetic first moment conversion, based on apparent conductivity depth map (CDI) andMeasure the reference first moment that magnetic field data obtains the abnormal area at acceptance point place;

Transient electromagnetic first moment transform definition is:

$M^{\left(1\right)}={\∫}_0^{\mathrm{\∞}}B\left(t\right)dt---\left(2\right)$

It is magnetic responsiveness was carved into ∞ integration from 0 o'clock. Resistive being limited in time domain of transient electromagnetic equals oneRank square, becomes resistive restricting data by magnetic-field measurement data transaction, is equal to first moment conversion.

Step C specifically comprises:

Sub-step C1:

Based on measurement data and even vertical magnetic field RESPONSE CALCULATION measurement data first moment greatly;

Sub-step C1 specifically comprises:

Because measurement data is in limited time range, in order to obtain the first moment of measurement data, needWant the magnetic field integral part outside polishing time range, here, the expression of the measurement data first moment of emulationFormula is:

Wherein, t₁And t_nRepresent respectively measurement data initial time and deadline, σ₁And σ_nTable respectivelyShow t₁And t_nThe apparent conductivity at place, head represents that magnetic field is from 0 to t₁Integration in time, middle part represents to surveyAmount data are from t₁To t_nIntegration, afterbody represents that magnetic field is from t_nArrive the integration in the ∞ time.

Middle part be the fitting function that first obtains measurement data B (t), by fitting function at time t₁And t_nInNumerical integration obtains middle part integration. Because the time range internal magnetic field at head and afterbody cannot be measured,Here adopt even magnetic responsiveness integration greatly in theory to replace, head is by calculating from 0 to t₁InEven magnetic responsiveness integration greatly obtains; Afterbody is by calculating from t_nRing to the even magnetic field greatly of ∞Answer integration to obtain, wherein even electrical conductivity is greatly estimated to obtain from apparent conductivity depth map.

In the measurement data first moment of above-mentioned emulation, relate to even magnetic responsiveness integration greatly in theorySpecifically solve as follows:

G (x, y, t) is the even vertical magnetic field response greatly of unitary current elementary excitation, and its expression formula is:

$\begin{array}{c}G\left(x,y,t\right)=I\frac{\mathrm{\μ}}{8{\mathrm{\πx}}^2}\[\frac{1}{{\mathrm{\γ}}^2}\left(2+\frac{x^2}{{\mathrm{\ρ}}^2}-2{\mathrm{\γ}}^2{x}^2\right)\frac{y}{\mathrm{\ρ}}erf\left(\mathrm{\γ}\mathrm{\ρ}\right)\\ {-\frac{2}{\sqrt{\mathrm{\π}}}\frac{x^2}{{\mathrm{\ρ}}^2}\frac{y}{\mathrm{\γ}}{e}^{\left(-{\mathrm{\γ}}^2{\mathrm{\ρ}}^2\right)}-\frac{2}{{\mathrm{\γ}}^2}{e}^{\left(-{\mathrm{\γ}}^2{x}^2\right)}erf\left(\mathrm{\γ}y\right)\]}_0^{\mathrm{\γ}\left(t\right)}\end{array}---\left(4\right)$

Wherein,σ represents even electrical conductivity greatly, from apparent conductivity depth map, estimatesObtain, μ is space permeability, and x and y represent the position of the end of the relative line current of corresponding acceptance point unitPut coordinate,T is receive path time window.

In the present embodiment, owing to adopting rectangle to determine greatly source loop line, therefore under rectangular loop source forcing,Even magnetic responsiveness B greatly_z(t) be:

$\begin{array}{c}{B}_z\left(t\right)={\[G(x,y,t)\]}_{(y_1,x_2,t)}^{(y_1,x_1,t)}+{\[G(x,y,t)\]}_{(y_2,x_1,t)}^{(y_2,x_2,t)}\\ +{\[G(x,y,t)\]}_{(x_1,y_2,,t)}^{(x_1,y_1,t)}+{\[G(x,y,t)\]}_{(x_2,y_2,t)}^{(x_2,y_1,t)}\end{array}---\left(5\right)$

In above formula, G (x, y, t) is the vertical magnetic field response of unitary current unit, x₁＝X_E-X_R，x₂＝X_W-X_R，y₁＝Y_N-Y_R，y₂＝Y_S-Y_R, wherein (X_E，Y_N) be the apex coordinate that rectangle is determined greatly source loop line northeastward,(X_W，Y_S) be the apex coordinate that rectangle is determined greatly source loop line southwestward, as shown in figure 12.

The integrated form of the vertical magnetic field response of definition unitary current unit:

$\begin{array}{c}K\left(x,y,t\right)={\∫}_t^{\mathrm{\∞}}G\left(x,y,\mathrm{\τ}\right)d\mathrm{\τ}\\ =I\frac{{\mathrm{\σ\μ}}^2}{64{\mathrm{\πx}}^3{\mathrm{\γ}}^4}\[\left(2+\frac{x^2}{{\mathrm{\ρ}}^2}-4{\mathrm{\γ}}^2{x}^2\right)\frac{y}{\mathrm{\ρ}}erf\left(\mathrm{\γ}\mathrm{\ρ}\right)\\ {-\frac{2}{\sqrt{\mathrm{\π}}}\frac{x^2}{{\mathrm{\ρ}}^2}{\mathrm{\γye}}^{\left(-{\mathrm{\γ}}^2{\mathrm{\ρ}}^2\right)}+2\left({\mathrm{\γ}}^2{x}^2-1\right)e^{\left(-{\mathrm{\γ}}^2{x}^2\right)}erf\left(\mathrm{\γ}y\right)\]}_0^{\mathrm{\γ}\left(t\right)}\\ +I\frac{{\mathrm{\σ\μ}}^2}{16\mathrm{\π}}x{\∫}_0^y\frac{erf\left(\mathrm{\γ}\mathrm{\ρ}\right)}{\mathrm{\ρ}}dy\end{array}---\left(6\right)$

In the time of t → 0, the expression formula of above formula bracket part equals zero, and the integral part of end can enterOne step abbreviation, therefore, rectangle is determined greatly the wherein table of the first moment of a frame excitation formation of source loop lineReaching formula is

$M_z^{\left(1\right)}(x,y_1,y_2)=I\frac{{\mathrm{\σ\μ}}^2}{16\mathrm{\π}}x\·{\[\ln (y+\mathrm{\ρ})\]}_{y_1}^{y_2}---\left(7\right)$

Above formula represents from frame (x, y₁) to (x, y₂) integration, determine greatly under the excitation of source loop line, all at rectangleThe first moment of even the earth equals the first moment sum of four edges frame.

Rectangle is determined greatly under frame excitation of source loop line, evenly large magnetic responsiveness from the time 0 tot₁Integration be:

${\∫}_0^{t_1}{B}_{{\mathrm{\σ}}_1}\left(t\right)dt={M_z}^{\left(1\right)}(x,y_1,y_2)-\[K(x,y_2,t_1)-K(x,y_1,t_1)\]---\left(8\right)$

Even magnetic responsiveness is greatly from time t_nIntegration to ∞ is:

${\∫}_{t_n}^{\mathrm{\∞}}{B}_{{\mathrm{\σ}}_n}\left(t\right)dt=K(x,y_2,t_n)-K(x,y_1,t_n)---\left(9\right)$

In the present embodiment, calculate rectangle according to expression formula (8) and determine greatly the even of every frame of source loop lineThe magnetic responsiveness of the earth is from the time 0 to t₁Integration, by the above-mentioned integrated value summation of four edges frame,Can try to achieve the head in measurement data first moment; Calculate rectangle according to expression formula (9) and determine greatly source loop lineThe even magnetic responsiveness greatly of every frame is from time t_nTo the integration of ∞, above-mentioned by four edges frameIntegrated value summation, can try to achieve the afterbody in measurement data first moment.

In order to make the measurement data of emulation be close to field survey data, by the measurement data of above-mentioned emulationFirst moment adds 5% white Gaussian noise;

Sub-step C2: estimated background electrical conductivity from apparent conductivity depth map (CDI), then utilizeEvenly the first moment expression formula of vertical magnetic field is obtained background response first moment greatly;

Sub-step C2 specifically comprises: from apparent conductivity depth map CDI, obtain initial time t₁With cutOnly time t_nThe apparent conductivity σ at place₁And σ_n, and estimate background conductance rate σ_bg。

Wherein above-mentioned estimation background conductance rate σ_bgSpecifically comprise:

Background response first moment can be divided and obtain from 0 to ∞ time inner product by even magnetic responsiveness greatly.

In the present embodiment, background response first moment, according to expression formula (7), is determined greatly source loop line by rectangleFour edges frame first moment be added obtain.

Sub-step C3: the ginseng of rejecting background first moment and obtain abnormal area from measurement data first momentExamine first moment;

Fig. 5 is that measurement data first moment, background first moment, abnormal area first moment are on each survey lineCurve map. The reference first moment of this abnormal area is as the reference first moment in ensuing invertingd＝(d₁，d₂，…，d_N)^T, in this example, receive the N=231 that counts.

Step D, the constraints of time constant vector in structure inverting; Constraint in described step DCondition is the one in following three kinds of constraintss:

Constraints a: in refutation process, the value in time constant vector be restricted on the occasion of, insteadDrill when initial, the initial value of time constant vector is obtained by apparent conductivity depth map;

Constraints b: in refutation process, the value in time constant vector is restricted on the occasion of, and rootApparent conductivity according to infinitesimal is given weighted value to the value in time constant vector, in the time that inverting is initial, timeBetween the initial value of constant vector be zero;

Constraints c: in refutation process, the value in time constant vector is restricted on the occasion of, and rootThe degree of depth according to infinitesimal is given weighted value to the value in time constant vector, when inverting is initial, and time constantThe initial value of vector is zero;

Wherein,

The initial value of the time constant vector in described constraints a is obtained by following mode:

First, definition initial time constant vectorThe time constant of k infinitesimalFor

${\mathrm{\τ}}_k^0=\frac{{\mathrm{\σ}}_k{\mathrm{\μL}}^2}{25.6}---\left(10\right)$

Wherein, σ_kBe the apparent conductivity value of k infinitesimal, by apparent conductivity depth map, CDI obtains,μ is space permeability, the length of side that L is infinitesimal.

According to least square method, the matching of actual measurement response and initial theoretical response is poor is

$S=\underset{n=1}{\overset{N}{\Σ}}{(d_n-\mathrm{\β}\·c_n^0)}^2---\left(11\right)$

Whereinβ is fitting coefficient, d_nFor the reference first moment of abnormal area. For makingThe poor minimum of matching, has partial differential equationSolve this partial differential equation, obtain β, when initialBetween constant vector be τ=β τ⁰。

In described constraints b, weighted value computational methods are:

The position in CDI is relevant to infinitesimal for electrical conductivity weighted value, and it is distributed between 0 and 1, theK infinitesimal electrical conductivity weighted value is

$W\left({\mathrm{\σ}}_k\right)=\frac{{\mathrm{\σ}}_k}{{\mathrm{\σ}}_{max}}---\left(12\right)$

Wherein, σ_kBe the apparent conductivity value of k infinitesimal, by apparent conductivity depth map, CDI obtains, σ_maxFor apparent conductivity maximum in whole CDI. Electrical conductivity weighting accounts for large electrical conductivity region in CDILeading role, weakens how much coupling factor impacts in little electrical conductivity region.

The weighted value computational methods of described constraints c are:

Depth weighted value is relevant to the degree of depth of infinitesimal, and the depth weighted value of k infinitesimal is

W(z_k，s₀，d₀)＝tanh(s₀·(z_k-d₀))(13)

In above formula, z_kBe the degree of depth of k infinitesimal, s₀For slope factor, d₀For reference depth. RightThe time constant of different depth infinitesimal is given weighted value, the depth weighted shadow that weakens shallow-layer coupling factorRing.

Step e, the reference first moment based on abnormal area and error thereof, the time constant that infinitesimal is correspondingThe object function of vector and how much coupling factor matrix structure complementary operations;

Step e specifically comprises:

c_nRepresent the theoretical first moment that n acceptance point place abnormal area is just being drilled computing, its computational methods are:

In theory, the total first moment response of the earth is folded by background response first moment and anomalous body response first momentAdd formation, expression formula is as follows:

$M_n^{\left(1\right)}={M_{0n}}^{\left(1\right)}+\underset{k=1}{\overset{K}{\Σ}}{G}_{nk}{\mathrm{\τ}}_k---\left(14\right)$

In above formula,Represent the background response first moment at n acceptance point place, the summation of equation the rightPart represents anomalous body response first moment, G_nkWhat be k infinitesimal with respect to n acceptance point is severalWhat coupling factor matrix, τ_kIt is the time constant of k infinitesimal. Total first moment of the earth in theory shouldEqual measurement data first moment,, need from measurement data first moment, be finally inversed by unknown τ here_k。

For simplified operation, in refutation process, background response first moment picks from measurement data first momentRemove, abnormal area is just being drilled the theoretical first moment of computing and is being got:

$c_n=\underset{k=1}{\overset{K}{\Σ}}{G}_{nk}{\mathrm{\τ}}_k---\left(15\right)$

The object function expression formula of complementary operation is

Wherein, the time constant of infinitesimal vector is τ=(τ₁，τ₂，…，τ_K)^T，τ_kIt is the time of k infinitesimalConstant, the error of the reference first moment of abnormal area is q=(q₁，q₂，…，q_N)^T，q_nBe n exceptions areaThe error of the reference first moment in territory, c_nRepresent that n acceptance point place abnormal area just drilling the theory one of computingRank square, G_nkBe how much coupling factor matrixes, d_nFor the reference first moment of abnormal area.

Step F, based on the constraints of time constant vector corresponding to infinitesimal, adopts optimization algorithmThe object function of complementary operation is carried out to iteration;

Step F specifically comprises:

Sub-step F1:

If adopt constraints a, the vector of the time constant before the i+1 time iteration is τⁱ, the i+1 timeTime constant vector after iteration is τⁱ⁺¹＝τⁱ+δτⁱ；

If adopt constraints b, the vector of the time constant before the i+1 time iteration is τⁱ, the i+1 timeTime constant vector after iteration is τⁱ⁺¹＝τⁱ+W(σ_k)δτⁱ；

If adopt constraints c, the vector of the time constant before the i+1 time iteration is τⁱ, the i+1 timeTime constant vector after iteration is τⁱ⁺¹＝τⁱ+W(z_k，s₀，d₀)δτⁱ；

Sub-step F2:

Adopt optimization algorithm to carry out iteration to described object function.

Optimization algorithm is the one in following methods: steepest descent method, conjugate gradient method, Newton method,Quasi-Newton method, Trust Region, singular value decomposition method, simulated annealing, genetic algorithm.

In the present embodiment, the explanation alternative manner as an example of steepest descent method example,

δτⁱCalculate by following formula:

${\mathrm{\δ\τ}}^i=-{\mathrm{\α}}^i\·\▿S\left({\mathrm{\τ}}^i\right)---\left(17\right)$

Step-length αⁱExpression formula be

${\mathrm{\α}}^i=-\frac{\underset{n=1}{\overset{N}{\Σ}}{e}_n^i{\mathrm{\ψ}}_n^i}{\underset{n=1}{\overset{N}{\Σ}}{\left({\mathrm{\ψ}}_n^i\right)}^2}---\left(18\right)$

In above formula, intermediate variableWithExpression formula be

$e_n^i=\frac{(d_n-c_n^i)}{q_n}---\left(19\right)$

${\mathrm{\ψ}}_n^i=\underset{k=1}{\overset{K}{\Σ}}\frac{1}{q_n}{G}_{nk}\·{\▿}_{k}S\left({\mathrm{\τ}}^i\right)---\left(20\right)$

Wherein,Be the gradient of the poor object function of matching, be expressed as

${\▿}_{k}S\left({\mathrm{\τ}}^i\right)=-\frac{2}{N}\underset{n=1}{\overset{N}{\Σ}}\left(\frac{d_n-c_n^i}{q_n}\right)\frac{1}{q_n}{G}_{nk}---\left(21\right)$

Step G, judges whether the object function after iteration restrains, if convergence, while preserving optimizationBetween constant vector as inverting final result, i.e. τ=(τ₁，τ₂，…，τ_K)^T, otherwise, return to step F and holdOK, wherein, τ_kThe corresponding time constant of electrical conductivity of k infinitesimal in abnormal area, according toThe inverting final result of time constant vector obtains the position of anomalous body.

Wherein, whether described object function restrains refers to whether object function is less than 1, is less than 1 receiptsHold back, otherwise do not restrain.

Timeconstantτ_kRelevant with the electrical conductivity of corresponding k infinitesimal, electrical conductivity is larger, and the time is normalNumber is larger.

Fig. 6 is the inverting time constant sectional drawing of unconfined condition, and this inversion result through 73 times repeatedlyGeneration, 60s consuming time. Due to unconfined condition in refutation process, there is negative value in inverting time constant, withTime, because how much coupling factors of shallow-layer are larger, abnormal area ground proximity in inversion result. VisibleNonrestraint inversion PSD can not reflect the electric structure truly of dull and stereotyped anomalous body, therefore, and in refutation process,Setting constraints is very important.

Fig. 7 is the inverting time constant of the initial model of CDI that relates to of the constraints a in step DSectional drawing.

The initial matching of the initial model of CDI is poor is 1.58, and iteration 8 times, reaches matching after time 0.5sPoor being less than is 1 convergence. As can be seen from Figure 7, time constant inversion result can convergence CDIForm, the time constant maximum of inverting and constant value actual time of dull and stereotyped anomalous body(τ_slab=4.5msec) very identical.

Fig. 8 is the sectional drawing of the electrical conductivity weighted value that relates to of the constraints b in step D.

Electrical conductivity weighting makes in CDI large electrical conductivity region account for leading role, weakens little electrical conductivity regionHow much coupling factor impacts.

Fig. 9 is the inverting time constant sectional drawing that adopts electrical conductivity weighting. Initial matching is poor is 67, anti-After drilling iteration 50 times, reach convergence, 43s consuming time, as shown in the figure, is different from the initial model of CDIInversion result, in electrical conductivity weighting inversion result, high time constant region is confined to the volume of actual plateIn scope, reflect more exactly position and the volume of anomalous body.

Figure 10 is the sectional drawing of the depth weighted value that relates to of the constraints c in step D.

In this example, s₀＝0.004，d₀=350, the depth weighted value that this setup parameter calculates is as figureShown in 10.

Figure 11 adopts depth weighted inverting time constant sectional drawing. Initial matching is poor is 67, repeatedlyAfter 5 times, reach convergence, 8s consuming time, as shown in the figure, high time constant region and dull and stereotyped anomalous bodyPosition and volume coincide preferably.

So far, by reference to the accompanying drawings the present embodiment be have been described in detail. Describe according to above, thisThose skilled in the art should be to the present invention is based on the transient electromagnetic quick three-dimensional inverting of various boundary conditionsMethod has had clearly understanding.

In addition, the above-mentioned definition to each element and method is not limited in the various tools of mentioning in embodimentBody structure, shape or mode, those of ordinary skill in the art can change simply or replace it.

In sum, the invention provides a kind of transient electromagnetic quick three-dimensional based on various boundary conditionsInversion method, the method has adopted a kind of processing method of data compression and the positive algorithm of the three-dimensional of simplification,Solve the data volume that current 3-d inversion faces and greatly, just drilled a complicated difficult problem; Because 3-d inversion is askedThat inscribes is less qualitative, and direct inversion generally can not provide rational result, in order to impel inversion result to press close toActual underground structure, this method has proposed different constraints, can comparatively accurately be finally inversed by undergroundAnomalous body. The method that the present invention proposes can realize quick three of subsurface anomaly target in common computerDimension inverting has important application prospect in real time instant electric discharge power transformation magnetic data interpretation work.

Above-described specific embodiment, carries out object of the present invention, technical scheme and beneficial effectFurther description, institute it should be understood that the foregoing is only specific embodiments of the invention and, be not limited to the present invention, within the spirit and principles in the present invention all, any repairing of doingProtection scope of the present invention changes, be equal to replacement, improvement etc., within all should be included in.

Claims (10)

1. the transient electromagnetic quick three-dimensional inversion method based on various boundary conditions, its feature existsIn, comprising:

Steps A: lay on the ground emitter and acceptance point, set subsurface anomaly region, by differentNormal region is divided into the infinitesimal of cube shaped, according to how much of emitter, acceptance point and infinitesimalHow much coupling factor matrixes of calculation of parameter;

Step B: emitter emission current signal, electric current closes has no progeny, and each acceptance point gathers magnetic fieldData, then adopt even Earth model, and the magnetic field data of collection is converted to apparent conductivity depth map;

Step C: according to the conversion of transient electromagnetic first moment, based on apparent conductivity depth map and measurement magnetic fieldThe reference first moment of the abnormal area at data acquisition acceptance point place;

Step D: the constraints that builds time constant vector in inverting;

Step e: the reference first moment based on abnormal area and error thereof, the time constant that infinitesimal is correspondingThe object function of vector and how much coupling factor matrix structure complementary operations;

Step F: based on the constraints of time constant vector corresponding to infinitesimal, adopt optimization algorithmThe object function of complementary operation is carried out to iteration; And

Step G: judge whether the object function after iteration restrains, if convergence, while preserving optimizationBetween constant vector as inverting final result, i.e. τ=(τ₁，τ₂，…，τ_k)^T, otherwise, return to step F and holdOK, wherein, τ_kThe time constant corresponding to electrical conductivity of k infinitesimal in abnormal area, according to timeBetween the inverting final result of constant vector obtain position and the volume of anomalous body.

2. transient electromagnetic quick three-dimensional inversion method according to claim 1, is characterized in that,Constraints in described step D is the one in following three kinds of constraintss:

Constraints a: in refutation process, the value in time constant vector be restricted on the occasion of, insteadDrill when initial, the initial value of time constant vector is obtained by apparent conductivity depth map;

Constraints b: in refutation process, the value in time constant vector is restricted on the occasion of, and rootApparent conductivity according to infinitesimal is given weighted value to the value in time constant vector, in the time that inverting is initial, timeBetween the initial value of constant vector be zero;

Constraints c: in refutation process, the value in time constant vector is restricted on the occasion of, and rootThe degree of depth according to infinitesimal is given weighted value to the value in time constant vector, when inverting is initial, and time constantThe initial value of vector is zero.

3. transient electromagnetic quick three-dimensional inversion method according to claim 2, is characterized in that,The initial value of the time constant vector in described constraints a is obtained by following mode

Definition initial time constant vectorThe time constant of k infinitesimal is:

${\mathrm{\τ}}_k^0=\frac{{\mathrm{\σ}}_k{\mathrm{\μL}}^2}{25.6}$

Wherein, σ_kBe the apparent conductivity value of k infinitesimal, obtained by apparent conductivity depth map, μ isSpace permeability, the length of side that L is infinitesimal;

According to least square method, the matching of actual measurement response and initial theoretical response is poor is:

$S=\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{(d_n-\mathrm{\β}\·c_n^0)}^2$

Whereinβ is fitting coefficient, d_nFor the reference first moment of abnormal area, G_nkForHow much coupling factor matrixes, for making the poor minimum of matching, have partial differential equationSolve this partially micro-Divide equation, obtain β, initial time constant vector is τ=β τ⁰。

4. transient electromagnetic quick three-dimensional inversion method according to claim 2, is characterized in that,In described constraints b, weighted value computational methods are:

Electrical conductivity weighted value is relevant to the position of infinitesimal in apparent conductivity depth map, and it is distributed in 0 HeBetween 1, k infinitesimal electrical conductivity weighted value is:

$W\left({\mathrm{\σ}}_k\right)=\frac{{\mathrm{\σ}}_k}{{\mathrm{\σ}}_{max}}$

Wherein, σ_kBe the apparent conductivity value of k infinitesimal, obtained σ by apparent conductivity depth map_maxForThe apparent conductivity of maximum in whole apparent conductivity depth map.

5. transient electromagnetic quick three-dimensional inversion method according to claim 2, is characterized in that,The weighted value computational methods of described constraints c are:

Depth weighted value is relevant to the degree of depth of infinitesimal, and the depth weighted value of k infinitesimal is:

W(z_k，s₀，d₀)＝tanh(s₀·(z_k-d₀))

In above formula, z_kBe the degree of depth of k infinitesimal, s₀For slope factor, d₀For reference depth.

6. transient electromagnetic quick three-dimensional inversion method according to claim 2, is characterized in that,Described step F specifically comprises:

Sub-step F1:

If adopt constraints a, the vector of the time constant before the i+1 time iteration is τⁱ, the i+1 timeTime constant vector after iteration is τⁱ⁺¹＝τⁱ+δτⁱ；

If adopt constraints b, the vector of the time constant before the i+1 time iteration is τⁱ, the i+1 timeTime constant vector after iteration is τⁱ⁺¹＝τⁱ+W(σ_k)δτⁱ；

If adopt constraints c, the vector of the time constant before the i+1 time iteration is τⁱ, the i+1 timeTime constant vector after iteration is τⁱ⁺¹＝τⁱ+W(z_k，s₀，d₀)δτⁱ; And

Sub-step F2:

Adopt optimization algorithm to carry out iteration to described object function.

7. according to the transient electromagnetic quick three-dimensional described in any one claim in claim 1 to 6Inversion method, is characterized in that, described steps A specifically comprises:

Sub-step A1: be detected region division emitter and acceptance point;

Sub-step A2: abnormal area is divided into infinitesimal, and wherein, this abnormal area is for being detected districtIn territory, may there is the target area of anomalous body; And

Sub-step A3: according to the geometric parameter computational geometry coupling of emitter, acceptance point and infinitesimalClose factor matrix.

8. transient electromagnetic quick three-dimensional inversion method according to claim 7, is characterized in that,How much coupling factor matrix G of described sub-step A3_nkFor:

$G_{nk}=B_{0,k}{V}_k\frac{3({\hat{b}}_k\×{\hat{r}}_{nk}){\hat{r}}_{nk}-{\hat{b}}_k}{4{{\mathrm{\πr}}_{nk}}^3}\frac{{\mathrm{\π}}^2}{10}$

In above formula,Represent that k infinitesimal center point to the unit direction vector of n acceptance point,Represent the unit direction vector of the primary field that incides k infinitesimal, V_kRepresent the body of k infinitesimalLong-pending, B_0，kRepresent to incide k infinitesimal primary field amplitude, rnk represents that k infinitesimal is to nThe distance of acceptance point.

9. transient electromagnetic quick three-dimensional inversion method according to claim 8, is characterized in that,Described step e specifically comprises:

The object function expression formula of complementary operation is:

$S\left(\mathrm{\τ}\right)=\frac{1}{N}\underset{n=1}{\overset{N}{\Σ}}{\left(\frac{d_n-c_n}{q_n}\right)}^2$

Wherein, the time constant of infinitesimal vector is τ=(τ₁，τ₂，…，τ_k)^T，τ_kIt is the time of k infinitesimalConstant, the error of the reference first moment of abnormal area is q=(q₁，q₂，…，q_N)^T，q_nBe n exceptions areaThe error of the reference first moment in territory, c_nRepresent that n acceptance point place abnormal area just drilling the theory one of computingRank square,G_nkBe how much coupling factor matrixes, d_nFor the reference first moment of abnormal area.

10. according to the transient electromagnetic quick three-dimensional described in any one claim in claim 1 to 6Inversion method, is characterized in that, described step C specifically comprises:

Sub-step C1: based on measurement data and even vertical magnetic field RESPONSE CALCULATION measurement data one greatlyRank square;

Sub-step C2: estimated background electrical conductivity from apparent conductivity depth map, then utilize evenly largeThe first moment expression formula of ground vertical magnetic field is obtained background response first moment; And

Sub-step C3: the ginseng of rejecting background first moment and obtain abnormal area from measurement data first momentExamine first moment.