Small loop source transient electromagnetic inversion method based on improved particle swarm optimization

Abstract

The invention discloses a small loop source transient electromagnetic inversion method based on an improved particle swarm algorithm, and belongs to the technical field of geophysical exploration. According to the method, the small loop transient electromagnetic detection device is adopted for underground detection, and on the basis of a traditional particle swarm algorithm, the position updating method of particles in each iteration is improved; different from the existing particle swarm algorithm applied to multilayer medium inversion, the method provided by the invention assumes the conductivity and depth of the kth layer of the current optimal estimation value as true values for calculation during iteration; the problem that the accuracy of a certain layer of parameters is lost due to the fact that the fitness of a whole objective function is reduced during iteration of a traditional particle swarm algorithm is effectively avoided. Compared with traditional smoke ring inversion, the method can effectively solve the problem that the dielectric layering boundary conductivity and layering height estimation are inaccurate; compared with a traditional particle swarm algorithm, the method is higher in noise immunity.

Description

Technical Field

[0001] The invention relates to the field of geophysical detection technology, and in particular to a small-loop source transient electromagnetic inversion method based on an improved particle swarm algorithm. Background Art

[0002] Transient electromagnetic method (TEM) is a geophysical detection method based on the principle of electromagnetic induction. It emits a pulse field underground through a loop or a long straight line source, and detects the decay characteristics of the induced secondary field of the underground medium over time after the primary pulse is turned off, thereby obtaining the physical and electrical parameters (conductivity, burial depth) of the underground medium.

[0003] Among them, small loop source detection devices are widely used in urban underground engineering detection, unexploded ordnance detection, coal mine aquifer detection and other fields due to their advantages of small size, light weight and relatively low cost. Therefore, in recent years, the processing and interpretation of small loop detection data has become a research focus in the field of transient electromagnetic detection.

[0004] The use of known physical and electrical parameters to calculate the attenuation characteristics of the secondary field is called forward modeling. Conversely, obtaining physical and electrical parameters through detection data and reconstructing underground structures is called inversion. Traditional TEM inversion mainly uses the "smoke ring" fast imaging method, which is an approximate processing method with fast imaging speed, but is insensitive to the depth of abnormal geological bodies and has low inversion accuracy. The transient electromagnetic inversion problem is inherently nonlinear. Therefore, in recent years, nonlinear global optimization algorithms such as particle swarm optimization (PSO), genetic algorithm, simulated annealing and artificial neural network have been widely used in the processing and interpretation of transient electromagnetic data, but their noise resistance needs to be further improved. Summary of the invention

[0005] In order to solve the problems of insensitivity to the depth of abnormal geological bodies and low inversion accuracy in the existing smoke ring inversion technology, the present invention proposes a small-loop source transient electromagnetic inversion method based on an improved particle swarm algorithm, which has stronger noise resistance than the traditional particle swarm algorithm.

[0006] To achieve the above object, the technical solution of the present invention is as follows:

[0007] A small loop source transient electromagnetic inversion method based on improved particle swarm algorithm, the method comprises the following steps:

[0008] Step 1, using a small loop transient electromagnetic detection device to emit a pulse on the surface to be measured, and using the secondary induced electromotive force received by the detection system as observation data;

[0009] Step 2: set the search range of underground medium parameters (such as conductivity, depth, etc.) according to the actual geological conditions, and regard the underground medium parameters as particle positions, and set the initial position of the particle population within the search range;

[0010] Step 3, the induced electromotive force data obtained by forward modeling of the particle position and the observation data are used to construct an objective function based on the L2 norm as the fitness of the particle swarm update, and then the particle position update of the particle swarm algorithm is executed based on the fitness of the particle; that is, the fitness of the particle is calculated based on the constructed new objective function, and then the improved particle swarm algorithm of the present invention is executed to update the particle position;

[0011] When the objective function value reaches the accuracy requirement or the number of iterations exceeds the expected setting, the current particle position can be considered to be the underground medium parameter obtained by transient electromagnetic inversion.

[0012] The method of the present invention is different from the traditional particle swarm algorithm and can avoid the problem of losing the accuracy of a certain layer of parameters in order to ensure the overall fitness reduction during the iteration process, thereby ensuring the inversion accuracy.

[0013] Furthermore, underground medium parameters include: conductivity and depth.

[0014] Furthermore, step 1 also includes: constructing an inverted initial layered geoelectric model, which includes several layers, each layer contains two underground medium parameters of conductivity and depth, and the depth of the innermost layer of the geoelectric model (medium layer depth) is defined as infinity.

[0015] Furthermore, in step 2, setting the initial position of the particle population within the search range includes:

[0016] Initialize the number of particles M, the position of each particle, the initial search speed when the particle is updated, the maximum number of iterations of the population maxcount, the maximum number of iterations of the segmented loop loopcount, the number of layers k of the geoelectric model that is not currently involved in the iteration, and the precision eps. The number of particles M means that there are M geoelectric models to be optimized, and the particle position represents the estimated value of the geoelectric parameter during optimization.

[0017] Furthermore, in step 3, the objective function constructed based on the L2 norm by combining the induced electromotive force data obtained by forward modeling of the particle position with the observed data is expressed as:

[0018]

[0019] Among them, V t is the induced voltage of the tth time channel in the actual measurement, V′ tis the theoretical value of the induced voltage of the tth time channel obtained by forward modeling with the estimated value of the parameter to be determined, and T is the total number of time channels during measurement. The objective function actually reflects the degree of fit between the measured value and the forward modeling value. The smaller the fitness value, the closer the forward modeling value is to the measured value, and the closer the current parameter estimate is to the true value of the underground medium parameter. Therefore, when the objective function is used as the fitness of the example individual, the historical minimum fitness value of the current particle individual corresponds to the individual's optimal position, and the historical minimum fitness value of the objective function of the particle group corresponds to the optimal position of the particle population.

[0020] Furthermore, in step 3, the particle position update of the particle swarm algorithm based on the fitness of the particles specifically includes:

[0021] Step 301, calculate the fitness value of each particle based on the objective function, and find the optimal position of the current particle and the particle group;

[0022] Step 302, the k-th layer parameters of the geoelectric model of the current optimal position of the group are regarded as theoretical values ​​and do not participate in the particle position update;

[0023] Step 303, calculating the inertia weight coefficient and the updated particle velocity;

[0024] Step 304, calculating the updated position of the particle based on the updated particle velocity, and then determining whether the individual particle position exceeds the search range, and if not, performing population migration;

[0025] Step 304, determine whether the particle fitness value meets the preset conditions, if so, update the optimal position of the particle.

[0026] Furthermore, in the calculation process of step 302, in order to improve the calculation accuracy, the original iteration method of the particle swarm algorithm is improved, and the overall cycle is improved into a segmented cycle. In each cycle, the currently obtained optimal value is used as a reference, and the parameters to be calculated contained in the kth layer are regarded as known true values, and do not participate in the update calculation of the particle position, while the parameters of the remaining layers are iteratively updated; after each subsequent cycle, the parameters to be calculated in the k+1th layer are regarded as known true values ​​in turn according to the optimal value obtained in the previous cycle, and the next calculation is performed, and it is iterated repeatedly until the accuracy condition is met or the maximum number of iterations is reached.

[0027] Furthermore, in step 303, the particle speed update formula is:

[0028] The velocity update formula of any mth particle at the iter+1 iteration is:

[0029]

[0030] in, Respectively represent the particle speed of the mth particle after iter+1 and iter iterations, r 1 、r 2 is a random number uniformly distributed in the interval (0,1), c 1 、c 2 is the learning factor. Preferably, the value is set to 2 during inversion. m is the optimal position of the individual, p g is the optimal position of the group, the initial value of the iteration number iter is a preset value, for example, set to 0, ω is the inertia weight, and its calculation formula is:

[0031] ω=0.99 iter rand() / 2+0.1

[0032] Among them, rand() is a random number uniformly distributed in the interval (0,1).

[0033] Furthermore, in step 304, the updated position of the particle is calculated based on the updated particle velocity as follows:

[0034] The position update formula of the mth particle at the iter+1 iteration is calculated by the following formula:

[0035]

[0036] in, They represent the particle positions of the mth particle after iter+1 and iter iterations respectively. r is the constraint factor. Preferably, r takes the value of 1 during inversion.

[0037] It should be noted that the particle position after iterative update needs to meet the position limit range. If it does not meet the limit, the position is reselected.

[0038] Furthermore, step 304 includes:

[0039] Calculate the fitness of the particles, determine whether the current fitness of each particle is less than the individual historical optimal value, if so, update the individual optimal value; determine whether the current fitness of each particle is less than the group historical optimal value, if so, update the group optimal value.

[0040] Furthermore, in step 3, when the objective function value reaches the accuracy requirement or the number of iterations exceeds the expected setting, it is considered that the current particle position is the underground medium parameter obtained by transient electromagnetic inversion, which specifically includes:

[0041] Determine whether the fitness corresponding to the current optimal value of the group meets the accuracy requirements or reaches the maximum number of iterations of the segmented loop. If not, continue the loop; if the requirements are met, save the results, get the optimal solution of the segmented loop, and determine whether the maximum number of iterations of the overall loop is reached; if not, increase k by 1, and then regard the parameters to be determined contained in the kth layer as known true values, do not participate in the update calculation of the particle position, and continue the next loop; if the requirements are met, save and output the results, and the obtained optimal position of the particle can be used as the true value of the inverted underground medium parameter.

[0042] The technical solution provided by the present invention brings at least the following beneficial effects:

[0043] Based on the traditional particle swarm algorithm, the present invention improves the method of updating the position of particles in each iteration. Different from the existing particle swarm algorithm applied to multi-layer medium inversion, the method proposed by the present invention assumes the conductivity and depth of the kth layer of the current optimal estimate as the true value for calculation during iteration, which effectively avoids the problem of the traditional particle swarm algorithm losing the accuracy of a certain layer of parameters due to the decrease in the fitness of the overall objective function during iteration. Compared with the traditional smoke ring inversion, the method of the present invention can effectively solve the problem of inaccurate estimation of the conductivity and layer height of the medium layer boundary; compared with the traditional particle swarm algorithm, the method of the present invention is more resistant to noise. BRIEF DESCRIPTION OF THE DRAWINGS

[0044] In order to more clearly illustrate the technical solutions in the embodiments of the present invention, the drawings required for use in the description of the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without creative work.

[0045] Figure 1 A flowchart of a small loop source transient electromagnetic inversion method based on an improved particle swarm algorithm provided in an embodiment of the present invention;

[0046] Figure 2 A schematic diagram of a three-layer geoelectric model used in an embodiment of the present invention;

[0047] Figure 3 It is a schematic diagram of the results of inversion under ideal conditions using the improved particle swarm algorithm, the particle swarm algorithm and the smoke ring algorithm respectively in the embodiment of the present invention;

[0048] Figure 4 It is a schematic diagram of the results of inversion using the improved particle swarm algorithm and the particle swarm algorithm respectively after adding random noise in the embodiment of the present invention. DETAILED DESCRIPTION

[0049] To make the purpose, technical solutions and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be described in detail and completely in conjunction with the drawings in the implementation of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present application, rather than all of the embodiments. Generally, the components of the embodiments of the present invention described and shown in the drawings can be arranged and designed using different configurations. Therefore, the following detailed description of the embodiments of the present invention provided in the drawings is not intended to limit the scope of the present application for protection, but merely represents selected embodiments of the present invention.

[0050] The embodiment of the present invention provides a small-loop source transient electromagnetic inversion method based on an improved particle swarm algorithm, which uses a small-loop transient electromagnetic detection device for underground detection. On the basis of the traditional particle swarm algorithm, the position update method of particles in each iteration is improved; different from the existing particle swarm algorithm applied to multi-layer medium inversion, the method of the embodiment of the present invention assumes the conductivity and depth of the kth layer of the current optimal estimate as the true value for calculation during iteration, effectively avoiding the problem of the traditional particle swarm algorithm losing the accuracy of a certain layer of parameters due to the decrease in the fitness of the overall objective function during iteration; compared with the traditional smoke ring inversion, the method of the embodiment can effectively solve the problem of inaccurate estimation of the conductivity and layer height of the medium layer boundary; compared with the traditional particle swarm algorithm, the method has stronger noise resistance.

[0051] As a possible implementation, Figure 1 As shown, the small loop source transient electromagnetic inversion method based on improved particle swarm algorithm provided by the embodiment of the present invention includes the following steps:

[0052] S1: Use a small loop transient electromagnetic detection device to emit a pulse on the surface to be measured, and receive the secondary induced electromotive force generated by the underground medium under the excitation of the pulse as observation data;

[0053] S2: construct the initial layered geoelectric model for inversion;

[0054] S3: Determine the particle search range and initialize the particle population based on actual geological conditions;

[0055] S4: Calculate the fitness value of each particle and find the optimal position of the current particle and the particle group;

[0056] S5: The k-th layer parameter of the current optimal position of the group is regarded as a theoretical value and does not participate in the particle position update;

[0057] S6: Calculate the inertia weight coefficient and the particle speed update formula;

[0058] S7: Determine whether the position of individual particles exceeds the search range, and the population migrates;

[0059] S8: Determine whether the particle fitness value meets the conditions and update the optimal position of the particle;

[0060] S9: Determine whether the iteration termination condition is met, if yes, terminate the inversion algorithm and output the result, otherwise return to S5.

[0061] Preferably, a three-layer geoelectric model can be used, such as Figure 2 As shown, each layer contains two parameters: conductivity and layer depth, σ 1 =σ 3 =0.01(S / m),σ 2 =0.05(S / m),h 1 =50m,h 2 =100m, the depth of the third dielectric layer is defined as infinite.

[0062] A small loop transient electromagnetic detection device is used to emit a pulse on the surface to be tested, and the secondary induced electromotive force V = (V 1 ,V 2 ,...,V T ) as the observation data, where T is the total number of time channels during measurement. In this embodiment, the detection device uses a coil with a radius of 0.75m and 50 turns, the emission current intensity is 10A, and the delay time is 10 -6 s to 10 0 s, with equal logarithmic intervals, a total of 30 time channels.

[0063] Set the number of particles M = 60, the maximum number of iterations of the population maxcount = 50, the maximum number of iterations of the segment loop loopcount = 10, the number of layers not currently involved in the iteration k = 0, and the precision eps = 1e-6. m =(σ 1 ,σ 2 ,σ 3 ,h 1 ,h 2 ) is defined as the particle position, the position range is limited to plus or minus 50% of the geoelectric parameter, and the initial particle group position is randomly selected within the position range. The initial search speed υ is a random number uniformly distributed in the (0,1) interval.

[0064] Calculate the induced electromotive force V' of the geoelectric model represented by each particle according to the forward formula m =(V′ 1 ,V′ 2 ,...,V′ T ), according to the objective function fitness formula Calculate the fitness value, and record the current particle position as the individual optimal position p mAmong the current M particles, the position of the particle with the smallest fitness is recorded as the optimal position p of the group. g .

[0065] The kth layer parameter of the current swarm optimal position is regarded as a theoretical value and does not participate in the particle position update (k = 0 means that all layer parameters need to be iteratively calculated, and it is reset to 1 when the k value exceeds the number of geoelectric layers).

[0066] According to the formula ω=0.99 iter rand() / 2+0.1 calculates the inertia weight coefficient.

[0067] According to the formula Calculates the update speed of the particle.

[0068] According to the formula Calculate the particle migration position and determine whether it exceeds the search range. If it exceeds the search range, re-value it.

[0069] Calculate the objective function fitness, determine whether the current objective function fitness of each particle is less than the individual historical optimal value, if so, update the individual optimal value; determine whether the current objective function fitness of each particle is less than the group historical optimal value, if so, update the group optimal value.

[0070] Determine whether the fitness of the objective function corresponding to the current optimal value of the group meets the accuracy requirements or reaches the maximum number of iterations of the segmented loop. If the conditions are not met, continue the loop; if the requirements are met, save the results, obtain the optimal solution of the segmented loop, and determine whether the maximum number of iterations of the overall loop is reached. If the conditions are not met, the parameters to be determined contained in the k=k+1 layer are regarded as known true values, and do not participate in the update calculation of the particle position, and continue to the next loop; if the requirements are met, save and output the results, and the obtained optimal position of the particle can be used as the true value of the inverted underground medium parameter.

[0071] The results of inversion using the improved particle swarm algorithm, the particle swarm algorithm and the smoke ring algorithm in the embodiments of the present invention are as follows: Figure 3 The error table is shown in Table 1.

[0072] Table 1 Comparison of inversion results (ideal case) (unit: σ (S / m), h (m))

[0073]

[0074] Depend on Figure 3 As shown in Table 1, the smoke ring inversion method can only invert the conductivity of the first layer, and the depth estimation is inaccurate, and the deep geoelectric parameters cannot be accurately obtained; while the inversion errors of the two nonlinear algorithms are both within 1%, and the inversion accuracy of the first layer of media is higher.

[0075] In the embodiment of the present invention, when the signal-to-noise ratio is 20 dB, the results of inversion using the improved particle swarm algorithm and the particle swarm algorithm are as follows: Figure 4 The error table is shown in Table 2.

[0076] Table 2 Comparison of inversion results (SNR = 20dB) (unit: σ (S / m), h (m))

[0077]

[0078] Depend on Figure 4 As can be seen from Table 2, after adding random noise, the overall inversion errors of the two algorithms have increased, which is mainly reflected in the geoelectric parameters of the second and third layers. Among them, the inversion error of the particle swarm algorithm for the conductivity of the third layer is 49.37%, which is close to the set boundary range of 50%, and the inversion error of the second layer depth also exceeds 30%. Compared with the particle swarm algorithm, the small loop source transient electromagnetic inversion method based on the improved particle swarm algorithm (referred to as improved PSO) proposed in the embodiment of the present invention controls all geoelectric parameter errors within 20% when adding random noise for inversion, and the inversion accuracy is significantly higher than that of the particle swarm algorithm ( Figure 4 The PSO in the algorithm has stronger noise immunity.

[0079] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, rather than to limit it. Although the present invention has been described in detail with reference to the aforementioned embodiments, those skilled in the art should understand that they can still modify the technical solutions described in the aforementioned embodiments, or make equivalent replacements for some of the technical features therein. However, these modifications or replacements do not deviate the essence of the corresponding technical solutions from the spirit and scope of the technical solutions of the embodiments of the present invention.

[0080] The above are only some embodiments of the present invention. For those skilled in the art, several modifications and improvements can be made without departing from the creative concept of the present invention, which all belong to the protection scope of the present invention.

Claims

1. A small loop source transient electromagnetic inversion method based on improved particle swarm algorithm, characterized in that: The following steps are involved: Step 1, using a small loop transient electromagnetic detection device to emit a pulse on the surface to be measured, and using the secondary induced electromotive force received by the detection system as observation data; Step 2: set the search range of underground medium parameters according to actual geological conditions, regard the underground medium parameters as particle positions, and set the initial position of the particle population within the search range; Step 3, the induced electromotive force data obtained by forward modeling of the particle position and the observed data are used to construct the objective function based on the L2 norm as the fitness of the particle swarm update, and then the particle position update of the particle swarm algorithm is executed based on the fitness of the particle; When the objective function value reaches the accuracy requirement or the number of iterations exceeds the expected setting, the current particle position is considered to be the underground medium parameter obtained by transient electromagnetic inversion.

2. The method according to claim 1, characterized in that The underground medium parameters include: conductivity and depth; and step 1 also includes: constructing an inversion initial layered geoelectric model, the geoelectric model includes several layers, each layer contains two underground medium parameters of conductivity and depth, and the depth of the innermost layer of the geoelectric model is defined as infinity.

3. The method according to claim 2, characterized in that In step 2, setting the initial position of the particle population within the search range includes: Initialize the number of particles M, the position of each particle, the initial search speed when updating particles, the maximum number of iterations of the population, the maximum number of iterations of the segmented loop, the number of layers k of the geoelectric model that is not currently involved in the iteration, and the accuracy.

4. The method according to claim 2, characterized in that In step 3, the objective function constructed based on the L2 norm by combining the induced electromotive force data obtained by forward modeling of the particle position with the observed data is expressed as: Among them, V t is the induced voltage of the tth time channel in the actual measurement, V′ t is the theoretical value of the induced voltage in the tth time channel obtained by forward modeling with the estimated values ​​of the parameters to be determined, and T is the total number of time channels during measurement.

5. The method according to any one of claims 3 to 4, characterized in that: In step 3, the particle position update of the particle swarm algorithm based on the fitness of the particles specifically includes: Step 301, calculate the fitness value of each particle based on the objective function, and find the optimal position of the current particle and the particle group; Step 302, the k-th layer parameters of the geoelectric model of the current optimal position of the group are regarded as theoretical values ​​and do not participate in the particle position update; Step 303, calculating the inertia weight coefficient and the updated particle velocity; Step 304, calculating the updated position of the particle based on the updated particle velocity, and then determining whether the individual particle position exceeds the search range, and if not, performing population migration; Step 304, determine whether the particle fitness value meets the preset conditions, if so, update the optimal position of the particle.

6. The method according to claim 2, characterized in that In step 302, the iteration mode of the particle swarm algorithm is improved from a whole cycle to a segmented cycle. In each segment of the cycle, the currently obtained optimal value is used as a reference, and the parameters to be determined contained in the kth layer are regarded as known true values ​​and do not participate in the update calculation of the particle position, while the parameters of the remaining layers are iteratively updated; After each cycle, the k+1th layer of the geoelectric model parameters to be calculated will be regarded as known true values ​​according to the optimal value obtained in the previous cycle, and the next calculation will be performed, and it will be iterated repeatedly until the accuracy condition is met or the maximum number of iterations is reached.

7. The method according to claim 2, characterized in that In step 303, the particle speed update formula is: in, They represent the particle speed of the mth particle after the iter+1th and iterth iterations, r1 and r2 are random numbers uniformly distributed in the (0,1) interval, c1 and c2 are two preset learning factors, and p m is the optimal position of the individual, p g is the optimal position of the group, the initial value of the iteration number iter is the preset value, ω is the inertia weight, and its calculation formula is: ω=0.99 iter ·rand() / 2+0.1 Among them, rand() is a random number uniformly distributed in the interval (0,1).

8. The method according to claim 6, characterized in that In step 304, the updated position of the particle is calculated based on the updated particle velocity as follows: in, They represent the particle positions of the mth particle after iter+1 and iter iterations respectively, and r is the preset constraint factor.

9. The method according to claim 2, characterized in that Step 304 includes: Calculate the fitness of the particles, determine whether the current fitness of each particle is less than the individual historical optimal value, if so, update the individual optimal value; determine whether the current fitness of each particle is less than the group historical optimal value, if so, update the group optimal value.

10. The method according to claim 2, characterized in that In step 3, when the objective function value reaches the accuracy requirement or the number of iterations exceeds the expected setting, it is considered that the current particle position is the underground medium parameter obtained by transient electromagnetic inversion, which specifically includes: Determine whether the fitness corresponding to the current optimal value of the group meets the accuracy requirements or reaches the maximum number of iterations of the segmented loop. If not, continue the loop; if satisfied, save the result, get the optimal solution of the segmented loop, and determine whether the maximum number of iterations of the overall loop is reached; if not satisfied, increase k by 1, and then regard the parameters to be determined contained in the kth layer as known true values, do not participate in the update calculation of the particle position, and continue the next loop; if satisfied, save and output the result, and the obtained optimal position of the particle can be used as the true value of the inverted underground medium parameter.

Images