CN114488314A - Geological inversion method based on land and underwater direct current combined measurement - Google Patents

Geological inversion method based on land and underwater direct current combined measurement Download PDF

Info

Publication number
CN114488314A
CN114488314A CN202111508309.8A CN202111508309A CN114488314A CN 114488314 A CN114488314 A CN 114488314A CN 202111508309 A CN202111508309 A CN 202111508309A CN 114488314 A CN114488314 A CN 114488314A
Authority
CN
China
Prior art keywords
electrode
water
inversion
data
land
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202111508309.8A
Other languages
Chinese (zh)
Other versions
CN114488314B (en
Inventor
王堃鹏
王向鹏
焦永春
李棂泽
纪伟
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Chengdu Univeristy of Technology
Original Assignee
Chengdu Univeristy of Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Chengdu Univeristy of Technology filed Critical Chengdu Univeristy of Technology
Priority to CN202111508309.8A priority Critical patent/CN114488314B/en
Publication of CN114488314A publication Critical patent/CN114488314A/en
Application granted granted Critical
Publication of CN114488314B publication Critical patent/CN114488314B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V3/00Electric or magnetic prospecting or detecting; Measuring magnetic field characteristics of the earth, e.g. declination, deviation
    • G01V3/02Electric or magnetic prospecting or detecting; Measuring magnetic field characteristics of the earth, e.g. declination, deviation operating with propagation of electric current
    • G01V3/04Electric or magnetic prospecting or detecting; Measuring magnetic field characteristics of the earth, e.g. declination, deviation operating with propagation of electric current using dc
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V3/00Electric or magnetic prospecting or detecting; Measuring magnetic field characteristics of the earth, e.g. declination, deviation
    • G01V3/38Processing data, e.g. for analysis, for interpretation, for correction
    • YGENERAL 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
    • Y02TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02ATECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/30Assessment of water resources

Landscapes

  • Life Sciences & Earth Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Environmental & Geological Engineering (AREA)
  • Geology (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

The invention discloses a geological inversion method based on land and underwater direct current combined measurement, which comprises the following steps: s1: arranging electrodes in a target region, wherein the target region comprises a water region, lands are respectively arranged on two sides of the water region, and the electrodes are arranged from the land on one side to the bottom of the water region step by step until the electrodes are arranged on the land on the other side; s2: determining the spatial distribution of the water body in the water area by adopting a sonar detection device, collecting water samples at different depths and different points in the water area for analysis, determining the resistivity spatial distribution structure of the water body, and establishing a three-dimensional model of the water body in the water area and the fluctuation of the water bottom topography; the invention is based on direct current combined measurement, and from the practical perspective, designs a plurality of technical key points, including: the geological condition from the shore to the deep water area is effectively controlled by a progressive measuring method from land to water; eliminating the influence of flowing water on the disturbance of the electrode; and establishing a reservoir model, and eliminating the influence of the water body on the inversion of the underground medium.

Description

Geological inversion method based on land and underwater direct current combined measurement
Technical Field
The invention belongs to the technical field of land and underwater direct current detection, and particularly relates to a geological inversion method based on land and underwater direct current combined measurement.
Background
So far, the forward and backward evolution method and theory of direct current are mature, and the ground direct current method has wide application and is relatively stable. However, in view of the current development, the methods for the combined measurement and inversion of terrestrial and underwater direct currents are rarely studied and the technology is still incomplete. The safety detection of reservoir and hydroelectric engineering dam bodies is limited by the limitation of a conventional method, and breakthrough is needed urgently, and the development of the research of the land and underwater direct current combined measurement and inversion method considering the dynamic water environment is particularly important.
Disclosure of Invention
The invention provides a geological inversion method based on land and underwater direct current combined measurement, and aims to solve the problems.
The invention is realized in such a way that a geological inversion method based on land and underwater direct current combined measurement comprises the following steps:
s1: arranging electrodes in a target region, wherein the target region comprises a water region, the two sides of the water region are respectively provided with lands, and the electrodes are arranged from the land on one side to the bottom of the water region step by step until the electrodes are arranged on the land on the other side;
s2: determining the spatial distribution of the water body in the water area by adopting a sonar detection device, collecting water samples of different depths and different points of the water area for analysis, determining the resistivity spatial distribution structure of the water body, and establishing a three-dimensional model of the water body in the water area and the topography fluctuation of the water bottom;
s3: according to the electrode arrangement on the land and the bottom of the water area, a temperature sodium device or a dipole-dipole device is adopted to change the positions of a power supply electrode and a receiving electrode, different power supply-receiving observation systems are established, and observation data are obtained;
s4: preprocessing the observation data, and removing abnormal data points in the observation data;
s5: acquiring an apparent resistivity pseudo-section diagram according to the preprocessed observation data, analyzing the apparent resistivity pseudo-section diagram, and preliminarily judging the underground geological condition in a target region;
s6: establishing different geological models according to the water body and geological characteristics in a target region, and performing two-dimensional and three-dimensional forward modeling simulation calculation analysis, wherein the whole forward modeling calculation region adopts an unstructured triangular unit or a tetrahedral unit to perform two-dimensional and three-dimensional mesh generation;
s7: according to the three-dimensional model of the water body and the water bottom topography fluctuation in the water area, an initial model containing a three-dimensional water body resistivity structure is established, and grid subdivision is carried out by using unstructured triangular units or tetrahedral units so as to realize two-dimensional and three-dimensional inversion.
Further, in step S3, the transforming the positions of the power supply electrode and the receiving electrode by using the wenner device or the dipole-dipole device specifically includes: adopting A, B, M, N four electrodes, wherein electrode A is anode, electrode B is cathode, and electrode M and electrode N are measuring point electrodes; continuously changing the positions of the electrode A, the electrode B, the electrode M and the electrode N, thereby inputting current to different depths underground and obtaining the potential difference of the electrode M, N at the current measuring point; recording the positions of the electrode A, the electrode B, the electrode M and the electrode N, the current and the voltage of each transformation, and finally calculating the apparent resistivity in the following mode:
Figure BDA0003405032650000021
wherein, Delta UMNM, N, measuring the potential difference between the electrode pair, I is the measuring current, K is the device coefficient, the specific calculation rules of the device coefficients K of the wenner device and the dipole-dipole device are respectively as follows:
Figure BDA0003405032650000022
Figure BDA0003405032650000023
wherein, AM, AN, BM, BN and MN are the distances between the electrodes respectively.
Further, in step S4, the preprocessing the observation data specifically includes: firstly, Fourier transform is carried out on observation data, disturbance of flowing water on the electrode position and influence of the flowing water on the data are analyzed, and filtering processing is carried out on the data of the same observation point along with time; and then carrying out Fourier inversion, superposing data of observation points changing along with time, solving the mean value and the variance, further sorting the data into apparent resistivity data, editing the apparent resistivity data of each point, and removing the data beyond the normal range.
Further, in step S5, the obtaining the apparent resistivity pseudo-cross-sectional view specifically includes: for the Wener device, an apparent resistivity pseudo-section diagram is obtained by taking a half of the distance between an electrode A and an electrode B as a longitudinal axis and taking a central point between an electrode M and an electrode N as a transverse position of a measuring point; for the dipole-dipole device, half of the distance between the middle point between the electrodes A, B and the middle point between the electrodes MN is the vertical axis, the central point of the measuring device is the transverse position of the measuring point, and the apparent resistivity profile is obtained.
Further, in step S6, the performing two-dimensional and three-dimensional forward modeling calculation analysis specifically includes: setting the uniform resistivity of a water body, the resistivity of an irregular abnormal body below the water bottom and the specific background resistivity of the rest parts;
the two-dimensional potential calculation equation and the three-dimensional potential calculation equation are respectively as follows:
-▽[σ(x,z)▽φ(x,y,z)]=I(x,y,z)
-▽[σ(x,y,z)▽φ(x,y,z)]=I(x,y,z)
wherein σ (x, z) is a two-dimensional distribution of the conductivity, σ (x, y, z) is a three-dimensional distribution of the conductivity, φ (x, y, z) is a three-dimensional distribution of the potential, and I (x, y, z) is a current source term;
using a finite element method to disperse the upper expression to obtain a linear equation system:
Kv=S
k is a sparse matrix and comprises grid unit geometric information and resistivity structure information, v is a vector to be solved, and S is a right-end item;
and further utilizing the vector v and the interpolation matrix to obtain the apparent resistivity value of any measuring point position in the calculation space.
Further, in step S7, the method specifically includes:
setting an inversion target function:
Figure BDA0003405032650000031
wherein m is a column vector composed of all resistivity parameters, λ is a regularization factor,
Figure BDA0003405032650000032
in order to be an objective function of the data,
Figure BDA0003405032650000046
is a model objective function;
Figure BDA0003405032650000041
where Δ d is the residual vector between the observation data and the forward data, WdWeighting the matrix for the data;
Wd=diag{1/σ1,1/σ2,...1/σj...1/σm}
wherein, the data weighting matrix is a diagonal matrix, and the diagonal element is the reciprocal of each data variance;
Figure BDA0003405032650000042
wherein, WmThe model covariance matrix is used to control the smoothness of the model in the transverse and longitudinal directions, and the matrix is used to constrain the local model using a gradient operator or a laplacian operator.
Further, in step S7, the inversion method adopts a gauss-newton method, which specifically includes:
s71: selecting an initial resistivity model m0Setting the maximum iteration times of inversion and the minimum fitting difference;
s72: solving a system of linear equations
Figure RE-GDA0003527008210000044
Obtaining an iterative update direction pk
S73: setting an initial step length, starting one-dimensional linear search, optimizing the model updating step length alpha of the current iteration, and updating the model according to the following formula:
mk+1=mk+αpk
s74: if the current fitting difference is smaller than the minimum threshold value or the iteration times are larger than the maximum iteration times, terminating inversion; if the conditions are not met, jumping to the second step and starting new iteration;
in the linear system of equations, the following sensitivity matrices are included:
Figure BDA0003405032650000044
the element of the sensitivity matrix is the partial derivative of different observation data about the model parameter, and the matrix is a full-rank matrix;
in the inversion procedure, the fitting difference is calculated as follows:
Figure BDA0003405032650000045
further, in the inversion process, λ is adjusted by the following mechanism:
when the inversion fitting difference is difficult to reduce, the lambda is reduced to one tenth of the original lambda;
when the linear search is difficult to find the appropriate update step length, λ is reduced to one tenth of the original.
Compared with the prior art, the invention has the beneficial effects that: the invention is based on direct current combined measurement, and from the practical perspective, designs a plurality of technical key points, including: the geological condition from the shore to the deep water area is effectively controlled by a progressive measuring method from land to water; eliminating the influence of flowing water on the disturbance of the electrode; and establishing a water reservoir model to eliminate the influence of the water body on the inversion of the underground medium.
Drawings
FIG. 1 is a schematic diagram of the general arrangement of land-underwater-land electrodes provided by the present invention;
FIG. 2 is a flow chart of an embodiment provided by the present invention;
FIG. 3 is a schematic diagram of an electrode arrangement provided by an embodiment of the present invention;
FIG. 4 is a schematic diagram of sonar detection provided by an embodiment of the present invention;
FIG. 5 is a schematic diagram of the movement of the power supply electrode and the measurement electrode according to the embodiment of the present invention;
FIG. 6 is a resistivity forward model provided by an embodiment of the present invention;
fig. 7 is a schematic diagram of forward mesh generation according to an embodiment of the present invention;
FIG. 8 is a graph of a cross-sectional apparent resistivity profile provided by an embodiment of the present invention;
FIG. 9 shows a first inversion result provided by an embodiment of the present invention;
FIG. 10 is an inverse mesh generation provided by an embodiment of the present invention;
FIG. 11 is a second inversion result provided in accordance with an embodiment of the present invention;
fig. 12 shows an inversion result three provided by the embodiment of the present invention.
Detailed Description
In order to make the objects, technical solutions and advantages of the present invention more apparent, the present invention will be described in further detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative of the invention and are not intended to limit the invention.
In the description of the present invention, it is to be understood that the terms "length", "width", "upper", "lower", "front", "rear", "left", "right", "vertical", "horizontal", "top", "bottom", "inner", "outer", etc. indicate orientations or positional relationships based on those shown in the drawings, and are merely for convenience of description and simplicity of description, but do not indicate or imply that the device or element referred to must have a particular orientation, be constructed in a particular orientation, and be operated, and thus, are not to be construed as limiting the present invention. In addition, in the description of the present invention, "a plurality" means two or more unless otherwise specifically limited.
Example 1
As shown in fig. 1, in order to achieve the purpose of gradually passing from the shore to the deep water, thereby effectively controlling the related geologic bodies under and around the water body, the invention firstly requires that the electrode arrangement is gradually arranged from the land to the water and is arranged on the land again;
in order to reduce the influence of a water body on inversion, the method requires that a sonar device is adopted before an observation system and inversion are designed, a three-dimensional model of the water body and the water bottom topography of the current whole work area is controlled, water samples at different depths and different point positions are collected, the three-dimensional distribution of the resistivity of the water body is analyzed, and whether obvious resistivity macroscopically layering exists or not is judged.
According to the land and underwater resistivity arrangement shown in fig. 1, the positions of the power supply electrode and the receiving electrode can be changed by adopting a wenner device or a dipole-dipole device and the like, and different power supply-receiving observation systems are established, so that more accurate geological recognition capability is achieved.
Due to the effects of current flow and subsurface currents on the data, abnormal data points must be removed from the observed data. Firstly, Fourier transform is carried out, the disturbance of moving water on the electrode position and the influence of data are analyzed, and filtering processing is carried out on the change data of the same observation point along with time. And then carrying out Fourier inversion, overlapping data of observation points changing along with time, solving the mean value and the variance, further sorting the data into apparent resistivity data, editing the apparent resistivity data of each point, and removing the data beyond the normal range.
For the Wener device, a horizontal position of a measuring point is taken as a center point between an electrode M and an electrode N, and an apparent resistivity pseudo-cross section diagram is obtained by taking a half of the distance between the electrode A and the electrode B as a longitudinal axis; for a dipole-dipole device, half of the distance between the middle point between the electrodes A, B and the middle point between the electrodes M, N is a vertical axis, the central point of the measuring device is the transverse position of a measuring point, and an apparent resistivity pseudo-section diagram is obtained. According to the imaging result, horizontal slices and vertical slices at different positions are cut out, and the position of the leakage point can be preliminarily judged. And an initial model is provided for forward simulation analysis and backward calculation, and inversion multi-solution is reduced.
In order to adapt to complex terrains, the invention requires the mesh subdivision to be carried out by using triangular or tetrahedral unstructured units. In addition, different geological models are established according to the water body and geological characteristics of the current work area, two-dimensional and three-dimensional forward modeling analysis is carried out, and the actual resolution capability of the current observation system is preliminarily judged.
In order to eliminate the influence of the water body on the inversion, the invention requires that an initial model containing a three-dimensional water body resistivity structure is established according to the water body and the water bottom three-dimensional structure, and a triangular or tetrahedral unstructured unit is used for mesh generation to realize two-dimensional and three-dimensional inversion. And marking the resistivity structure and the topographic relief of the water body before inversion, and taking the resistivity structure and the topographic relief as fixed parameters in the inversion without participating in model modification.
Based on direct current combined measurement, the invention designs a plurality of technical key points from the practical angle, which comprises the following steps: the geological condition from the shore to the deep water area is effectively controlled by a progressive measuring method from land to water; eliminating the influence of the flowing water on the disturbance of the electrode; and establishing a reservoir model, and eliminating the influence of the water body on the inversion of the underground medium.
Example 2
As shown in fig. 2, the invention discloses a geological inversion method based on land and underwater direct current combined measurement, which comprises the following steps:
1. land-underwater-land electrode arrangement
As shown in FIG. 3, the vertical direction Z represents elevation and the horizontal direction X represents horizontal distance. The electrodes are required to be arranged from the land to the water and then to the land, the underwater electrodes are in contact with the water bottom surface, the water depth is about 1.5 m, and the electrodes pass through the water from the left bank and are arranged to the right bank.
2. Establishing water model
As shown in fig. 4, according to sonar detection results, water body space distribution is established; acquiring the actual water surface and water body elevation according to GPS measurement; and determining the resistivity spatial distribution structure of the water body assembly according to water sample analysis of different point positions and depths of the water area.
If the water body has obvious macroscopic resistivity stratification, a water body resistivity structure with a stratified structure needs to be established so as to remove the influence of water in the forward inversion more accurately.
3. Establishing an observation system
According to the land and underwater resistivity arrangement shown in fig. 1, the positions of the power supply electrode and the receiving electrode can be changed by adopting a wenner device, a dipole-dipole device and the like, and different power supply-receiving observation systems are established, so that more accurate geological recognition capability is achieved.
As shown in fig. 5, the winner device is applied to geological exploration, and the usage of the winner device is to continuously change the positions of the anode a, the cathode B and the measuring point electrode M, N, so as to input current to different depths in the underground and obtain the current potential difference of MN. The A, B, M, N position, current and voltage of each transformation need to be recorded, and the final apparent resistivity is calculated as follows:
Figure BDA0003405032650000081
wherein, Delta UMNTo measure M, N the potential difference between the electrode pair, I the measured current, and K the device coefficient, the specific calculation rule is:
Figure BDA0003405032650000082
wherein, AM, AN, MN are the distance between every electrode respectively.
4. Data pre-processing
And removing the forward data with higher sensitivity to stabilize the inversion calculation of the next step. If the processed data is actual measurement data, abnormal flying spots in the actual measurement data need to be analyzed, and the actual measurement voltage and current need to be superposed. Through spectrum analysis and filtering, the influence of flowing water and other extraneous currents on the observed data can be removed.
In addition, through multiple observations, the variance of the current data point needs to be calculated, and the data with larger variance has smaller weight in the inversion, so that the effect of the data with poor quality in the inversion is reduced.
5. Qualitative analysis
The electrode arrangement mode is a temperature-sensitive device, half of the distance between the anode A (+) and the cathode B (-) is taken as a longitudinal axis, the central point between the measuring electrode pair M and the measuring electrode pair N is taken as the transverse position of a measuring point, and an apparent resistivity section drawing is obtained, wherein the apparent resistivity is calculated by taking the formula as the basis.
By analyzing the apparent resistivity profile, the underground geological conditions, especially the lateral positions of potential leakage points, can be preliminarily judged.
6. Simulation calculation of grid section and forward modeling
The forward modeling is an important component in the inversion calculation and is also a key analysis process of the practical problem. Therefore, mesh generation and forward modeling calculation need to be realized in step 6.
Fig. 6 is a forward modeling diagram provided in the embodiment, in which the uniform resistivity of the water body is set to be 10 Ω m, the resistivity of the irregular anomaly below the water bottom is set to be 20 Ω m, the resistivity of the background of the rest part is set to be 100 Ω m, and the electrode distribution is consistent with that in fig. 3.
As shown in fig. 7, the whole forward calculation area is subdivided by using unstructured triangular units, and the triangular units have strong capability of adapting to complex terrain, so that the calculation accuracy of forward inversion can be improved.
The potential calculation equation is as follows:
-▽[σ(x,z)▽φ(x,y,z)]=I(x,y,z)
wherein σ (x, z) is the two-dimensional distribution of the conductivity, ∑ φ (x, y, z) is the three-dimensional distribution of the potential, and I (x, y, z) is a current source term;
using a finite element method to disperse the upper expression to obtain a linear equation system:
Kv=S
k is a sparse matrix and comprises grid unit geometric information and resistivity structure information, v is a vector to be solved, and S is a right-end item;
the linear equation set can be solved by direct solution or iterative solution, the direct solution precision is high, and when the problem to be solved is small, the direct solution is fast. The iterative solution is suitable for large-scale calculation and has low memory consumption.
And further utilizing the final solution vector v and the interpolation matrix to obtain the apparent resistivity value of any measured point position in the calculation space.
Fig. 8 is a final apparent resistivity cross-section of the model of fig. 6, from which it can be seen that the distribution of apparent resistivity may reflect to some extent the presence of anomalies.
7. Inverse calculation
The whole calculation area is subdivided by adopting an unstructured triangulation unit, the water body and land structures are discretized, and the resistivity values of different model areas, the total number of electrodes, the geometric distance between the electrodes, the electrode water penetration depth, the type of an electrode device, the total number of apparent resistivities, the model size and the maximum inversion iteration number can be set.
The inversion calculation of the invention firstly needs to set a target function, and the inversion target function in the embodiment is as follows:
Figure BDA0003405032650000101
wherein m is a column vector composed of all resistivity parameters, λ is a regularization factor,
Figure BDA0003405032650000107
in order to be an objective function of the data,
Figure BDA0003405032650000108
is a model objective function;
Figure BDA0003405032650000102
wherein, Δ d is a residual vector of the observation data and the forward data, the observation data and the forward data are logarithmic visual resistivity values, WdWeighting the matrix for the data;
Wd=diag{1/σ1,1/σ2,...1/σj...1/σm}
wherein, the data weighting matrix is a diagonal matrix, and the diagonal element is the reciprocal of each data variance; the method has the effects of controlling the weight of different data in inversion, improving the influence of high-quality data and reducing low-quality data. In addition, the data weighting matrix can also play a role in normalization, and is very helpful for stabilizing inversion.
Figure BDA0003405032650000103
Wherein, WmThe model covariance matrix is used to control the smoothness of the model in the transverse and longitudinal directions, and the matrix is used to constrain the local model using a gradient operator or a laplacian operator.
λ is a regularization factor, which acts as a control
Figure BDA0003405032650000104
And
Figure BDA0003405032650000105
the relationship between them. The larger the lambda value is, the smoother the reverse resistivity structure is, and the lower the resolution is; the smaller the lambda is, the higher the resolution of the inversion result is, but the action of the data item in the inversion is stronger than that of the model item, so that the multi-solution property of the integral inversion is enhanced, and a redundant error structure appears in the inversion. For the inversion of the measured data, the lambda can be selected from empirical values of 0.01 to 10, and the reliability of the inversion result can be verified through artificial setting before the inversion.
In the inversion process, λ is adjusted by the following mechanism:
when the inversion fitting difference is difficult to reduce, the lambda is reduced to one tenth of the original lambda;
when the linear search is difficult to find the appropriate update step length, λ is reduced to one tenth of the original.
Before the inversion starts, an initial model needs to be set, and the closer to the initial model of the real model, the more accurate the inversion can be obtained. If the inversion data is actually measured data, resistivity values of the water body and the surrounding rock need to be obtained according to well logging, geological drilling and the like, so that the resistivity distribution of the known structure is controlled, and the requirement of accurate inversion is met.
The inversion method selected by the invention is specifically a Gauss-Newton method, and the flow of the Gauss-Newton method is as follows: :
s71: selecting an initial resistivity model m0Setting the maximum iteration times of inversion and the minimum fitting difference;
s72: solving a system of linear equations
Figure BDA0003405032650000106
Obtaining an iterative update direction pk
S73: setting an initial step length, starting one-dimensional linear search, optimizing the model updating step length alpha of the current iteration, and updating the model according to the following formula:
mk+1=mk+αpk
s74: if the current fitting difference is smaller than the minimum threshold value or the iteration times are larger than the maximum iteration times, terminating inversion; if the conditions are not met, jumping to the second step and starting new iteration;
it should be noted that step optimization degenerates the current inversion into a univariate optimization problem. Further, in the linear equation system, the following sensitivity matrix is included:
Figure BDA0003405032650000111
the element of the sensitivity matrix is the partial derivative of different observation data about the model parameter, and the matrix is a full-rank matrix; when the solution problem is large, the matrix is slow to solve. Therefore, the invention uses the reciprocity theorem to solve the sensitivity matrix so as to accelerate the inversion calculation.
In the inversion procedure, the fitting difference is calculated as follows:
Figure BDA0003405032650000112
based on the inversion process and steps, 2% of random noise is added on the basis of forward data of the model shown in FIG. 6 to serve as observation data, and the inversion effect of the method is verified.
Firstly, the water body is used as a free parameter for inversion, the water body and land structures are inverted from 100 Ω m, the inversion result is shown in fig. 9, and the inversion grid subdivision is shown in fig. 10. As can be seen from the results of fig. 9, even in underwater exploration, if the water body is inverted as a free parameter and the initial resistivity is consistent with that of the land, it is difficult to detect the low-resistance anomaly below the water bottom. This is because the attraction of the current to the electrokinetic field is very strong for bodies of water versus underground anomalies. Therefore, the water body is used as a free parameter to participate in inversion, underwater abnormity can be suppressed, and the purpose of exploring the underwater geological structure cannot be completed.
In order to avoid the above problems, the water body is used as a known structure to participate in the inversion, that is, the water body structure distribution is obtained according to the sonar detection result, but the water body is still used as a free parameter and is not fixed in the inversion, and the final inversion result is shown in fig. 11. From the results in the figures, it can be seen that although the water body is already a known value, it is still difficult to identify the low-resistance anomaly at the water bottom as a self-parameter participating in the calculation.
In order to further improve the inversion effect, the third inversion calculation is carried out, and the water body is used as a known structure and is used as a fixed parameter in the whole process without modification. Fig. 12 shows the result of the inversion again, from which it can be seen that the water body is taken as a known and fixed parameter to participate in the inversion, which has a significant improvement in the low resistance recognition below the water bottom, and the inversion result shows that there is a large anomaly in the transverse 10m to 20m position, which substantially corresponds to the position of the underwater anomaly shown in fig. 6.
The present invention is not limited to the above preferred embodiments, but includes all modifications, equivalents, and improvements within the spirit and scope of the present invention.

Claims (8)

1. A geological inversion method based on land and underwater direct current combined measurement is characterized by comprising the following steps:
s1: arranging electrodes in a target region, wherein the target region comprises a water region, lands are respectively arranged on two sides of the water region, and the electrodes are arranged from the land on one side to the bottom of the water region step by step until the electrodes are arranged on the land on the other side;
s2: determining the spatial distribution of the water body in the water area by adopting a sonar detection device, collecting water samples at different depths and different points in the water area for analysis, determining the resistivity spatial distribution structure of the water body, and establishing a three-dimensional model of the water body in the water area and the fluctuation of the water bottom topography;
s3: according to the electrode arrangement on the land and the bottom of the water area, a temperature sodium device or a dipole-dipole device is adopted to change the positions of a power supply electrode and a receiving electrode, different power supply-receiving observation systems are established, and observation data are obtained;
s4: preprocessing the observation data, and removing abnormal data points in the observation data;
s5: acquiring an apparent resistivity pseudo-cross-section diagram according to the preprocessed observation data, analyzing the apparent resistivity pseudo-cross-section diagram, and preliminarily judging the underground geological condition in a target region;
s6: establishing different geological models according to the water body and geological characteristics in a target region, and performing two-dimensional and three-dimensional forward modeling simulation calculation analysis, wherein the whole forward modeling calculation region adopts unstructured triangular units or tetrahedral units to perform two-dimensional and three-dimensional mesh generation;
s7: according to the three-dimensional model of the water body and the water bottom topography fluctuation in the water area, an initial model containing a three-dimensional water body resistivity structure is established, and grid subdivision is carried out by using unstructured triangular units or tetrahedral units so as to realize two-dimensional and three-dimensional inversion.
2. The geological inversion method based on land and underwater direct current combined measurement as claimed in claim 1, wherein in step S3, the transforming the positions of the power supply electrode and the receiving electrode by using a wenner device or a dipole-dipole device specifically comprises: adopting A, B, M, N four electrodes, wherein electrode A is positive electrode, electrode B is negative electrode, and electrode M and electrode N are measuring point electrodes; continuously changing the positions of the electrode A, the electrode B, the electrode M and the electrode N, thereby inputting current to different depths underground and obtaining the potential difference of the electrode M, N at the current measuring point; recording the positions of the electrode A, the electrode B, the electrode M and the electrode N, the current and the voltage of each transformation, and finally calculating the apparent resistivity in the following mode:
Figure RE-FDA0003527008200000021
wherein, Delta UMNM, N, measuring the potential difference between the electrode pair, I is the measuring current, K is the device coefficient, the specific calculation rules of the device coefficients K of the wenner device and the dipole-dipole device are respectively as follows:
Figure RE-FDA0003527008200000022
Figure RE-FDA0003527008200000023
wherein, AM, AN, BM, BN and MN are the distances between the electrodes respectively.
3. The geological inversion method based on land and underwater direct current combined measurement as claimed in claim 1, wherein in step S4, the preprocessing of the observation data specifically comprises: firstly, Fourier transform is carried out on observation data, the disturbance of moving water on the electrode position and the influence of data are analyzed, and filtering processing is carried out on the change data of the same observation point along with time; and then carrying out Fourier inversion, superposing data of observation points changing along with time, solving the mean value and the variance, further sorting the data into apparent resistivity data, editing the apparent resistivity data of each point, and removing the data beyond the normal range.
4. The geological inversion method based on land and underwater direct current combined measurement as claimed in claim 2, wherein in step S5, the obtaining of the apparent resistivity pseudo-section map specifically comprises: for the Wener device, an apparent resistivity pseudo-section diagram is obtained by taking a half of the distance between an electrode A and an electrode B as a longitudinal axis and taking a central point between an electrode M and an electrode N as a transverse position of a measuring point; for the dipole-dipole device, half of the distance between the middle point between the electrodes A, B and the middle point between the electrodes M, N is a vertical axis, the central point of the measuring device is the transverse position of the measuring point, and an apparent resistivity pseudo-section diagram is obtained.
5. The geological inversion method based on the combined measurement of land and underwater direct current according to claim 1, wherein the performing of the two-dimensional and three-dimensional forward modeling calculation analysis in step S6 specifically comprises: setting the uniform or laminar resistivity of a water body, the resistivity of irregular abnormal bodies below the water bottom and the specific background resistivity of the rest parts;
the two-dimensional potential calculation equation and the three-dimensional potential calculation equation are respectively as follows:
Figure RE-FDA0003527008200000031
Figure RE-FDA0003527008200000032
wherein σ (x, z) is a two-dimensional distribution of the electrical conductivity, σ (x, y, z) is a three-dimensional distribution of the electrical conductivity,
Figure RE-FDA0003527008200000033
i (x, y, z) is a current source term;
using a finite element method to disperse the upper expression to obtain a linear equation system:
Kv=S
k is a sparse matrix and comprises grid unit geometric information and resistivity structure information, v is a vector to be solved, and S is a right-end item;
and further utilizing the vector v and the interpolation matrix to obtain the apparent resistivity value of any measuring point position in the calculation space.
6. The geological inversion method based on the combined measurement of land and underwater direct current according to claim 1, wherein in step S7, the method specifically comprises:
setting an inversion target function:
Figure RE-FDA0003527008200000034
wherein m is a column vector composed of all resistivity parameters, λ is a regularization factor,
Figure RE-FDA0003527008200000035
in order to be an objective function of the data,
Figure RE-FDA0003527008200000036
is a model objective function;
Figure RE-FDA0003527008200000037
where Δ d is the residual vector between the observation data and the forward data, WdWeighting the matrix for the data;
Wd=diag{1/σ1,1/σ2,...1/σj...1/σm}
wherein, the data weighting matrix is a diagonal matrix, and the diagonal element is the reciprocal of each data variance;
Figure RE-FDA0003527008200000038
wherein, WmThe model covariance matrix is used to control the smoothness of the model in both the horizontal and vertical directions, and the matrix uses a gradient operator or laplacian operator to constrain the local model.
7. The geological inversion method based on the land and underwater direct current combined measurement as claimed in claim 6, wherein in step S7, the inversion method adopts the gauss-newton method, which specifically comprises:
s71: selecting an initial resistivity model m0Setting the maximum iteration times of inversion and the minimum fitting difference;
s72: solving a system of linear equations
Figure RE-FDA0003527008200000041
Obtaining an iterative update direction pk
S73: setting an initial step length, starting one-dimensional linear search, optimizing the model updating step length alpha of the current iteration, and updating the model according to the following formula:
mk+1=mk+αpk
s74: if the current fitting difference is smaller than the minimum threshold value or the iteration times are larger than the maximum iteration times, terminating inversion; if the conditions are not met, jumping to the second step and starting new iteration;
in the linear system of equations, the following sensitivity matrices are included:
Figure RE-FDA0003527008200000042
the element of the sensitivity matrix is the partial derivative of different observation data about the model parameter, and the matrix is a full-rank matrix;
in the inversion procedure, the fitting difference is calculated as follows:
Figure RE-FDA0003527008200000043
8. a geological inversion method based on land and underwater direct current combined measurement according to claim 7, characterized in that in the inversion process, λ is adjusted by the following mechanism:
when the inversion fitting difference is difficult to reduce, the lambda is reduced to one tenth of the original lambda;
when the linear search is difficult to find the appropriate update step length, λ is reduced to one tenth of the original.
CN202111508309.8A 2021-12-10 2021-12-10 Geological inversion method based on land and underwater direct current combined measurement Active CN114488314B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111508309.8A CN114488314B (en) 2021-12-10 2021-12-10 Geological inversion method based on land and underwater direct current combined measurement

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111508309.8A CN114488314B (en) 2021-12-10 2021-12-10 Geological inversion method based on land and underwater direct current combined measurement

Publications (2)

Publication Number Publication Date
CN114488314A true CN114488314A (en) 2022-05-13
CN114488314B CN114488314B (en) 2022-12-23

Family

ID=81492717

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111508309.8A Active CN114488314B (en) 2021-12-10 2021-12-10 Geological inversion method based on land and underwater direct current combined measurement

Country Status (1)

Country Link
CN (1) CN114488314B (en)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117218258A (en) * 2023-11-08 2023-12-12 山东大学 Shield geological and tunnel visualization method, system, medium and equipment

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20130197891A1 (en) * 2012-01-31 2013-08-01 Michael L. Jessop Subsurface hydrogeologic system modeling
AU2014309129A1 (en) * 2013-08-20 2016-03-10 Technoimaging, Llc Systems and methods for remote electromagnetic exploration for mineral and energy resources using stationary long-range transmitters
CN106199742A (en) * 2016-06-29 2016-12-07 吉林大学 A kind of Frequency-domain AEM 2.5 ties up band landform inversion method
CN107632322A (en) * 2017-08-01 2018-01-26 安徽理工大学 A kind of cable system and exploitation method suitable for waters electrical prospecting
CN109461359A (en) * 2018-11-16 2019-03-12 高军 A kind of aqueous geological structure forward probe method in tunnel
US20210094660A1 (en) * 2018-01-24 2021-04-01 Ocean Floor Geophysics Inc. Devices, methods, and systems for underwater surveying

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20130197891A1 (en) * 2012-01-31 2013-08-01 Michael L. Jessop Subsurface hydrogeologic system modeling
AU2014309129A1 (en) * 2013-08-20 2016-03-10 Technoimaging, Llc Systems and methods for remote electromagnetic exploration for mineral and energy resources using stationary long-range transmitters
CN106199742A (en) * 2016-06-29 2016-12-07 吉林大学 A kind of Frequency-domain AEM 2.5 ties up band landform inversion method
CN107632322A (en) * 2017-08-01 2018-01-26 安徽理工大学 A kind of cable system and exploitation method suitable for waters electrical prospecting
US20210094660A1 (en) * 2018-01-24 2021-04-01 Ocean Floor Geophysics Inc. Devices, methods, and systems for underwater surveying
CN109461359A (en) * 2018-11-16 2019-03-12 高军 A kind of aqueous geological structure forward probe method in tunnel

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
SHUNGUO WANG ET AL.: "Joint inversion of lake-floor electrical resistivity tomography and boat-towed radio-magnetotelluric data illustrated on synthetic data and an application from the Aspo Hard Rock Laboratory site, Sweden", 《GEOPHYSICAL JOURNAL INTERNATIONAL》 *
吴小平等: "基于非结构网格的电阻率三维带地形反演", 《地球物理学报》 *
王堃鹏: "张量CSAMT三维主轴各向异性正反演研究", 《中国优秀博硕士学位论文全文数据库(博士) 基础科学辑》 *
蒋甫伟: "一种走航式快速探测海床土电阻率的技术方法研究", 《中国优秀博硕士学位论文全文数据库(硕士) 基础科学辑》 *

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117218258A (en) * 2023-11-08 2023-12-12 山东大学 Shield geological and tunnel visualization method, system, medium and equipment
CN117218258B (en) * 2023-11-08 2024-03-22 山东大学 Shield geological and tunnel visualization method, system, medium and equipment

Also Published As

Publication number Publication date
CN114488314B (en) 2022-12-23

Similar Documents

Publication Publication Date Title
WO2022227206A1 (en) Fully convolutional neural network-based magnetotelluric inversion method
CN113221393B (en) Three-dimensional magnetotelluric anisotropy inversion method based on non-structural finite element method
CN112949134B (en) Earth-well transient electromagnetic inversion method based on non-structural finite element method
CN105549106B (en) A kind of gravity multiple solutions inversion method
Zhou Electrical resistivity tomography: a subsurface-imaging technique
CN109209338B (en) Electrical observation system and detection method for detecting abnormities beside well
US11782183B2 (en) Magnetotelluric inversion method based on fully convolutional neural network
CN114896564B (en) Transient electromagnetic two-dimensional Bayesian inversion method adopting adaptive Thiessen polygon parameterization
Jiang et al. Nonlinear inversion of electrical resistivity imaging using pruning Bayesian neural networks
CN114488314B (en) Geological inversion method based on land and underwater direct current combined measurement
CN115238550A (en) Self-adaptive unstructured grid landslide rainfall geoelectric field numerical simulation calculation method
CN116187168A (en) Method for improving water depth inversion accuracy based on neural network-gravity information wavelet decomposition
Friedel Estimation and scaling of hydrostratigraphic units: application of unsupervised machine learning and multivariate statistical techniques to hydrogeophysical data
Başokur et al. Object-based model verification by a genetic algorithm approach: Application in archeological targets
Qiang et al. Optimized arrays for electrical resistivity tomography survey using Bayesian experimental design
CN111723491B (en) Grounding parameter acquisition method based on arbitrary layered soil green function
Demirel et al. Two‐dimensional joint inversions of cross‐hole resistivity data and resolution analysis of combined arrays
CN117092702A (en) Construction method of hole-tunnel induced polarization water detection structure and inversion water detection method
Zhou et al. Stochastic structure-constrained image-guided inversion of geophysical data
CN117933049A (en) Aviation transient electromagnetic three-dimensional inversion method based on deep learning
Feng et al. Contrast between 2D inversion and 3D inversion based on 2D high-density resistivity data
CN115563791B (en) Magnetotelluric data inversion method based on compressed sensing reconstruction
CN114114438B (en) Quasi-three-dimensional inversion method for ground-to-air transient electromagnetic data of loop source
Yang et al. Wireline logs constraint borehole-to-surface resistivity inversion method and water injection monitoring analysis
Tsourlos et al. Efficient 2D inversion of long ERT sections

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