CN113742987B - Direct-current resistivity joint inversion method based on least square-particle swarm strategy - Google Patents
Direct-current resistivity joint inversion method based on least square-particle swarm strategy Download PDFInfo
- Publication number
- CN113742987B CN113742987B CN202111054380.3A CN202111054380A CN113742987B CN 113742987 B CN113742987 B CN 113742987B CN 202111054380 A CN202111054380 A CN 202111054380A CN 113742987 B CN113742987 B CN 113742987B
- Authority
- CN
- China
- Prior art keywords
- resistivity
- inversion
- model
- objective function
- obs
- 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
Links
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/25—Design optimisation, verification or simulation using particle-based methods
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2111/00—Details relating to CAD techniques
- G06F2111/06—Multi-objective optimisation, e.g. Pareto optimisation using simulated annealing [SA], ant colony algorithms or genetic algorithms [GA]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/02—Reliability analysis or reliability optimisation; Failure analysis, e.g. worst case scenario performance, failure mode and effects analysis [FMEA]
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
- Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
- Y02A90/30—Assessment of water resources
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
The invention discloses a direct-current resistivity joint inversion method based on a least square-particle swarm strategy, which adopts a direct-current resistivity method to carry out geological exploration on a selected area so as to obtain actual observation data d obs The method comprises the steps of carrying out a first treatment on the surface of the First, an initial resistivity model m is constructed 0 Resistivity forward calculation is performed by using a finite element method to obtain forward response f 1 (m); forward response f to be obtained 1 (m) and actual observed data d obs Substituting the inversion resistivity data set into an initial inversion objective function, and inverting the resistivity of the direct current method by using a least square method to obtain an inversion resistivity data set m 1 The method comprises the steps of carrying out a first treatment on the surface of the Will invert the resistivity dataset m 1 As a new initial resistivity model m 2 Solving the updated inversion objective function phi by using a particle swarm algorithm 2 (m); finally, adopting the updated inversion objective function phi 2 And (m) carrying out inversion detection by a direct current resistivity method to obtain the accurate position of the cavity, thereby obtaining the size, the extending direction and the burial state of the karst geological development zone.
Description
Technical Field
The invention relates to the technical field of electrical geological exploration, in particular to a direct-current resistivity joint inversion method based on a least square-particle swarm strategy.
Background
Karst geological development zone is always a relatively complex and difficult-to-handle disaster geological problem affecting engineering construction foundation construction and building stability, and in engineering construction, particularly large-scale and important engineering, the existence of karst can generate a plurality of special bad building conditions and engineering stability problems, so that the position and the form of underground cavities are very necessary to be ascertained. However, due to complex karst morphology and irregular interfaces, the detection is greatly disturbed, and finally the detection difficulty is high.
At present, there are many karst cave detection modes, wherein the direct current resistivity method is widely applied to cave body detection due to strong anti-interference capability, high efficiency and abundant information content of the direct current resistivity method, but the interpretation of detection results of the direct current resistivity method becomes more difficult due to the non-uniqueness of solutions existing in geophysical exploration, and finally, the accuracy and resolution of karst landforms obtained by inversion are low, and the positions of the cavities in the karst landforms cannot be accurately obtained. Therefore, how to provide a method can effectively improve the accuracy and resolution of inversion interpretation, so that the position of a cavity can be accurately obtained; meanwhile, the size, the extending direction and the burial state of the karst geological development zone can be ascertained, and the karst geological characteristics can be determined, so that the method is a research direction in the industry.
Disclosure of Invention
Aiming at the problems in the prior art, the invention provides a direct current resistivity joint inversion method based on a least square-particle swarm strategy, which can effectively improve the accuracy and resolution of inversion interpretation, thereby accurately obtaining the position of a cavity; meanwhile, the size, the extending direction and the burial state of the karst geological development zone can be ascertained, and the karst geological characteristics can be determined.
In order to achieve the above purpose, the technical scheme adopted by the invention is as follows: a direct current resistivity joint inversion method based on a least square-particle swarm strategy comprises the following specific steps:
firstly, performing geological exploration on a selected area by adopting a direct-current resistivity method to obtain actual observation data d obs ；
Step two, firstly, constructing an initial resistivity model m 0 Resistivity forward calculation is performed by using a finite element method to obtain forward response f 1 (m)；
Step three, forward response f obtained in the step two is obtained 1 (m) and the actual observed data d obtained in step one obs Substituting the inversion resistivity data set into an initial inversion objective function, and inverting the resistivity of the direct current method by using a least square method to obtain an inversion resistivity data set m 1 The method comprises the steps of carrying out a first treatment on the surface of the The specific initial inversion objective function is:
φ 1 (m)＝α 1 ||f 1 (m)-d obs ||+β 1 ||m 1 -m 0 ||
wherein phi is 1 (m) is the inversion objective function, α 1 、β 1 F is the damping coefficient 1 (m) is a forward response, d obs For actual observation data, ||m 1 -m 0 The I is a model constraint term;
step four, inverting the resistivity data set m 1 As a new initial resistivity model m 2 Solving the updated inversion objective function phi by using a particle swarm algorithm 2 (m)；
Step five, finally adopting the updated inversion objective function phi obtained in the step four 2 And (m) carrying out inversion detection by a direct current resistivity method to obtain the accurate position of the cavity, thereby obtaining the size, the extending direction and the burial state of the karst geological development zone.
Further, the step of performing resistivity forward calculation by adopting a finite element method in the step two is as follows:
construction of functional equations
Wherein: n (N) i A form function; u (u) i Is the potential value on the triangle node;
converting the electric field boundary problem into a transformation problem to solve;
wherein: sigma is the conductivity value of the triangular unit; Γ is the integration boundary; a is a point source; omega is the selected area; n is the external normal direction of the boundary; u is the potential; r is the distance from the boundary to the origin;
the simplified functional expressions of the simultaneous formulas (1) and (2) are as follows:
wherein: k is a node stiffness matrix;
using matrix p= (0, …, I A ，…，0) T And performing variation on the formula (3) to obtain a linear equation set:
KU＝P (4)
solving the equation set of the formula (4) to obtain the potential u of each node, namely, forward response f (m).
Further, in the fourth step, the updated inversion objective function phi is solved by using a particle swarm algorithm 2 The specific steps of (m) are as follows:
a) Setting model search space [ m ] min ,m max ]The number of particles is L, the maximum iteration number is N, and the minimum fitting error is delta;
b) For new initial resistivity model m 2 Initializing its model parameters M and velocity parameters v, and calculating forward response f of the model 2 (m)；
c) Forward response f in step b) 2 (m) and the actual observed data d obtained in step one obs Substituting the model into the following inversion objective function, calculating the inversion objective function, and determining an individual optimal model pbest and a group optimal model gbest
φ 2 (m)＝α 2 ||f 2 (m)-d obs ||+β 2 ||m 2 -m 0 || (6)
Wherein phi is 2 (m) is the inversion objective function, α 2 、β 2 F is the damping coefficient 2 (m) is the forward response of the model, d obs For actual observation data, ||m 2 -m 0 The I is a constraint item of the current model;
d) Updating the individual optimal model parameters M and the speed parameters v according to the individual optimal model pbest and the group optimal model gbest determined in the step c);
e) According to the individual optimal model parameters M and the speed parameters v obtained in the step c), calculating forward model response again, calculating an inversion objective function, and determining an individual optimal model pbest and a group optimal model gbest again;
f) If the fit error of the group optimal model gbest is smaller than delta or the iteration number reaches N, terminating the iteration, and outputting an updated inversion objective function phi 2 (m); otherwise, returning to the step d) to continue iteration until the iteration termination condition is met.
Compared with the prior art, the method adopts the combination of the least square and particle swarm strategy to carry out direct current resistivity inversion, firstly establishes a resistivity model, and uses a nonlinear least square algorithm to carry out inversion during primary inversion; on the basis, a new resistivity model is constructed, inversion is carried out by adopting a particle swarm algorithm, and the advantages of simplicity in implementation, high searching efficiency and rapid convergence are combined to obtain an inversion optimal solution. The invention suppresses the problem of multiple solutions when solving inversion equations by using an electric sounding method, and improves the interpretation precision and resolution of DC sounding exploration data. Inversion problems of the direct current method are multi-parameter, multi-extremum and non-linear, and the essence is to find the optimal solution in the functional space of the observed data and the forward model parameters. Therefore, the invention can effectively improve the accuracy and resolution of inversion interpretation, thereby accurately obtaining the position of the cavity; meanwhile, the size, the extending direction and the burial state of the karst geological development zone can be ascertained, and the karst geological characteristics can be determined, so that reliable geological data basis is provided for solving the potential hazards of karst geology disasters and designing engineering construction schemes.
Drawings
FIG. 1 is an overall flow chart of the present invention;
FIG. 2 is a graph of inversion results from a DC resistivity method using least squares;
FIG. 3 is a graph of inversion results of a DC resistivity method using particle swarms;
FIG. 4 is a graph of the direct current resistivity method joint inversion results based on the least squares-particle swarm strategy, using the present invention.
Detailed Description
The present invention will be further described below.
As shown in fig. 1, the specific steps of the present invention are:
firstly, performing geological exploration on a selected area by adopting a direct-current resistivity method to obtain actual observation data d obs ；
Step two, firstly, constructing an initial resistivity model m 0 Resistivity forward calculation is performed by using a finite element method to obtain forward response f 1 (m); the method for forward calculation of the resistivity by adopting the finite element method comprises the following steps:
construction of functional equations
Wherein: n (N) i A form function; u (u) i Is the potential value on the triangle node;
converting the electric field boundary problem into a transformation problem to solve;
wherein: sigma is the conductivity value of the triangular unit; Γ is the integration boundary; a is a point source; omega is the selected area; n is the external normal direction of the boundary; u is the potential; r is the distance from the boundary to the origin;
the simplified functional expressions of the simultaneous formulas (1) and (2) are as follows:
wherein: k is a node stiffness matrix;
using matrix p= (0, …, I A ，…，0) T And performing variation on the formula (3) to obtain a linear equation set:
KU＝P (4)
solving the equation set of the formula (4) so as to obtain the potential u of each node, namely forward response f (m);
step three, forward response f obtained in the step two is obtained 1 (m) and the actual observed data d obtained in step one obs Substituting the inversion resistivity data set into an initial inversion objective function, and inverting the resistivity of the direct current method by using a least square method to obtain an inversion resistivity data set m 1 The method comprises the steps of carrying out a first treatment on the surface of the The specific initial inversion objective function is:
φ 1 (m)＝α 1 ||f 1 (m)-d obs ||+β 1 ||m 1 -m 0 || (5)
wherein phi is 1 (m) is the inversion objective function, α 1 、β 1 F is the damping coefficient 1 (m) is a forward response, d obs For actual observation data, ||m 1 -m 0 The I is a model constraint term;
step four, inverting the resistivity data set m 1 As a new initial resistivity model m 2 Solving the updated inversion objective function phi by using a particle swarm algorithm 2 (m); wherein the updated inversion objective function phi is solved by using a particle swarm algorithm 2 The specific steps of (m) are as follows:
a) Setting model search space [ m ] min ,m max ]The number of particles is L, the maximum iteration number is N, and the minimum fitting error is delta;
b) For new initial resistivity model m 2 Initializing its model parameters M and velocity parameters v, and calculating forward response f of the model 2 (m)；
c) Forward response f in step b) 2 (m) and the actual observed data d obtained in step one obs Substituting the model into the following inversion objective function, calculating the inversion objective function, and determining an individual optimal model pbest and a group optimal model gbest
φ 2 (m)＝α 2 ||f 2 (m)-d obs ||+β 2 ||m 2 -m 0 || (6)
Wherein phi is 2 (m) is the inversion objective function, α 2 、β 2 F is the damping coefficient 2 (m) is the forward response of the model, d obs For actual observation data, ||m 2 -m 0 The I is a constraint item of the current model;
d) Updating the individual optimal model parameters M and the speed parameters v according to the individual optimal model pbest and the group optimal model gbest determined in the step c);
e) According to the individual optimal model parameters M and the speed parameters v obtained in the step c), calculating forward model response again, calculating an inversion objective function, and determining an individual optimal model pbest and a group optimal model gbest again;
f) If the fit error of the group optimal model gbest is smaller than delta or the iteration number reaches N, terminating the iteration, and outputting an updated inversion objective function phi 2 (m); otherwise, returning to the step d) to continue iteration until the iteration termination condition is met;
step five, finally adopting the updated inversion objective function phi obtained in the step four 2 And (m) carrying out inversion detection by a direct current resistivity method to obtain the accurate position of the cavity, thereby obtaining the size, the extending direction and the burial state of the karst geological development zone.
As can be seen from fig. 2 to 4, the inversion result obtained by the method can effectively improve the accuracy and resolution of inversion interpretation, thereby accurately obtaining the position of the cavity; meanwhile, the size, the extending direction and the burial state of the karst geological development zone can be ascertained, and the karst geological characteristics can be determined.
The foregoing is only a preferred embodiment of the invention, it being noted that: it will be apparent to those skilled in the art that various modifications and adaptations can be made without departing from the principles of the present invention, and such modifications and adaptations are intended to be comprehended within the scope of the invention.
Claims (1)
1. The direct-current resistivity joint inversion method based on the least square-particle swarm strategy is characterized by comprising the following specific steps of:
firstly, performing geological exploration on a selected area by adopting a direct-current resistivity method to obtain actual observation data d obs ；
Step two, firstly constructing initial resistivityModel m 0 Resistivity forward calculation is performed by using a finite element method to obtain forward response f 1 (m) specifically:
construction of functional equations
Wherein: n (N) i A form function; u (u) i Is the potential value on the triangle node;
converting the electric field boundary problem into a transformation problem to solve;
wherein: sigma is the conductivity value of the triangular unit; Γ is the integration boundary; a is a point source; omega is the selected area; n is the external normal direction of the boundary; u is the potential; r is the distance from the boundary to the origin;
the simplified functional expressions of the simultaneous formulas (1) and (2) are as follows:
wherein: k is a node stiffness matrix;
using matrix p= (0,) i. A ，...，0) T And performing variation on the formula (3) to obtain a linear equation set:
KU＝P (4)
solving the equation set of the formula (4) so as to obtain the potential u of each node, namely forward response f (m);
step three, forward response f obtained in the step two is obtained 1 (m) and the actual observed data d obtained in step one obs Substituting the inversion resistivity data set into an initial inversion objective function, and inverting the resistivity of the direct current method by using a least square method to obtain an inversion resistivity data set m 1 The method comprises the steps of carrying out a first treatment on the surface of the The specific initial inversion objective function is:
φ 1 (m)＝α 1 ||f 1 (m)-d obs ||+β 1 ||m 1 -m 0 ||
wherein phi is 1 (m) is the inversion objective function, α 1 、β 1 F is the damping coefficient 1 (m) is a forward response, d obs For actual observation data, ||m 1 -m 0 The I is a model constraint term;
step four, inverting the resistivity data set m 1 As a new initial resistivity model m 2 Solving the updated inversion objective function phi by using a particle swarm algorithm 2 (m) specifically:
a) Setting model search space [ m ] min ,m max ]The number of particles is L, the maximum iteration number is N, and the minimum fitting error is delta;
b) For new initial resistivity model m 2 Initializing its model parameters M and velocity parameters v, and calculating forward response f of the model 2 (m)；
c) Forward response f in step b) 2 (m) and the actual observed data d obtained in step one obs Substituting the model into the following inversion objective function, calculating the inversion objective function, and determining an individual optimal model pbest and a group optimal model gbest
φ 2 (m)＝α 2 ||f 2 (m)-d obs ||+β 2 ||m 2 -m 0 ||
Wherein phi is 2 (m) is the inversion objective function, α 2 、β 2 F is the damping coefficient 2 (m) is the forward response of the model, d obs For actual observation data, ||m 2 -m 0 The I is a constraint item of the current model;
d) Updating the individual optimal model parameters M and the speed parameters v according to the individual optimal model pbest and the group optimal model gbest determined in the step c);
e) According to the individual optimal model parameters M and the speed parameters v obtained in the step c), calculating forward model response again, calculating an inversion objective function, and determining an individual optimal model pbest and a group optimal model gbest again;
f) If the fit error of the group optimal model gbest is smaller than delta or the iteration number reaches N, terminating the iteration, and outputting an updated inversion objective function phi 2 (m); otherwise, returning to the step d) to continue iteration until the iteration termination condition is met;
step five, finally adopting the updated inversion objective function phi obtained in the step four 2 And (m) carrying out inversion detection by a direct current resistivity method to obtain the accurate position of the cavity, thereby obtaining the size, the extending direction and the burial state of the karst geological development zone.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111054380.3A CN113742987B (en) | 2021-09-09 | 2021-09-09 | Direct-current resistivity joint inversion method based on least square-particle swarm strategy |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111054380.3A CN113742987B (en) | 2021-09-09 | 2021-09-09 | Direct-current resistivity joint inversion method based on least square-particle swarm strategy |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN113742987A CN113742987A (en) | 2021-12-03 |
| CN113742987B true CN113742987B (en) | 2023-06-27 |
Family
ID=78737415
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202111054380.3A Active CN113742987B (en) | 2021-09-09 | 2021-09-09 | Direct-current resistivity joint inversion method based on least square-particle swarm strategy |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN113742987B (en) |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106526329A (en) * | 2016-11-10 | 2017-03-22 | 广东电网有限责任公司电力科学研究院 | Method and device for measuring resistivity from land surface to deep earth |
| CN108204944A (en) * | 2018-01-13 | 2018-06-26 | 福州大学 | The Buried Pipeline rate prediction method of LSSVM based on APSO optimizations |
| CN109597136A (en) * | 2018-11-27 | 2019-04-09 | 中煤科工集团西安研究院有限公司 | A kind of mine total space transient electromagnetic data processing method |
-
2021
- 2021-09-09 CN CN202111054380.3A patent/CN113742987B/en active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106526329A (en) * | 2016-11-10 | 2017-03-22 | 广东电网有限责任公司电力科学研究院 | Method and device for measuring resistivity from land surface to deep earth |
| CN108204944A (en) * | 2018-01-13 | 2018-06-26 | 福州大学 | The Buried Pipeline rate prediction method of LSSVM based on APSO optimizations |
| CN109597136A (en) * | 2018-11-27 | 2019-04-09 | 中煤科工集团西安研究院有限公司 | A kind of mine total space transient electromagnetic data processing method |
Non-Patent Citations (4)
| Title |
|---|
| Inversion for magnetotelluric data using the particle swarm optimization and regularized least squares;Yi-an Cui等;Journal of Applied Geophysics;第181卷;第1-7页 * |
| 基于最小二乘的矿井电法超前探测联合反演方法研究;李飞等;煤矿安全;第41-44页 * |
| 基于粒子群的全空间瞬变电磁二维反演方法研究;黄炜伟;工矿自动化;第47卷(第4期);第79-84页 * |
| 海洋可控源电磁与地震一维联合储层参数反演;徐凯军;杜润林;刘展;;石油地球物理勘探;第51卷(第01期);第197-203页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN113742987A (en) | 2021-12-03 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Lan et al. | A high‐order fast‐sweeping scheme for calculating first‐arrival travel times with an irregular surface | |
| Zhuang et al. | Fracture modeling using meshless methods and level sets in 3D: framework and modeling | |
| Chen et al. | A method of DEM construction and related error analysis | |
| CN111638556B (en) | Magnetotelluric forward modeling method and device based on geospatial solution strategy and storage medium | |
| CN113051779B (en) | Numerical simulation method of three-dimensional direct-current resistivity method | |
| CN112116708B (en) | Method and system for obtaining three-dimensional geological entity model | |
| CN111948712B (en) | Pre-stack linear inversion method based on depth domain seismic record | |
| CN112949134A (en) | Earth-well transient electromagnetic inversion method based on non-structural finite element method | |
| CN108108579B (en) | Boundary processing method of coupling finite element method in direct-current resistivity element-free method | |
| You et al. | An Eulerian method for computing the coherent ergodic partition of continuous dynamical systems | |
| CN112163332A (en) | 2.5-dimensional natural electric field numerical simulation method based on infinite element coupling of natural units | |
| CN114970290A (en) | Ground transient electromagnetic method parallel inversion method and system | |
| Budd et al. | The geometric integration of scale-invariant ordinary and partial differential equations | |
| CN113742987B (en) | Direct-current resistivity joint inversion method based on least square-particle swarm strategy | |
| Dong et al. | Improvement of affine iterative closest point algorithm for partial registration | |
| Martinho et al. | On reducing the slope parameter in terrain-following numerical ocean models | |
| Adcroft | Numerical algorithms for use in a dynamical model of the ocean | |
| CN117148460A (en) | Efficient inversion method of grounded long-line source transient electromagnetic method | |
| CN115906559B (en) | Geoelectromagnetic adaptive finite element forward modeling method based on mixed grid | |
| Chen | A digital-discrete method for smooth-continuous data reconstruction | |
| CN113204905B (en) | Numerical simulation method for finite element by contact induced polarization method | |
| Lv et al. | A Kriging interpolation-based boundary face method for 3D potential problems | |
| Abas et al. | Finite element method for fluid structure interaction with hp-adaptivity | |
| You et al. | VIALS: An Eulerian tool based on total variation and the level set method for studying dynamical systems | |
| Zhao et al. | New method for estimating strike and dip based on structural expansion orientation for 3D geological modeling |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |