Underground conductor three-dimensional imaging method based on magnetoresistivity method inversion technology

Abstract

The invention discloses an underground conductor three-dimensional imaging method based on a magnetoresistivity method inversion technology. The method comprises the following steps: acquiring magnetic field three-direction component data of each observation point on an observation surface on site; preprocessing the three-direction component data of the magnetic field, and extracting the three-direction component data of the magnetic field for inversion; establishing a magnetic field three-direction component joint inversion data fitting difference function; establishing a multi-source information constraint inversion model objective function; fitting a difference function and a multi-source information constraint inversion model objective function through magnetic field three-direction component joint inversion data, and constructing a Tikhonov regularization inversion objective function; the current density distribution value corresponding to the minimum value of the Tikhonov regularization inversion objective function is solved, a current density map of the geologic body space is drawn, and the area with the maximum current density is the space position of the underground conductor. The stability and reliability of an inversion imaging result are improved, and the method is applied to positioning and detection of underground conductors.

Description

Technical Field

[0001] The present invention relates to the field of geophysical exploration technology, and in particular to a three-dimensional imaging method for underground conductors based on magnetoresistivity inversion technology. Background Art

[0002] The magnetoresistivity method is an emerging geophysical exploration method that artificially supplies a non-inductive (DC or low-frequency AC) current between two points, measures the magnetic field it generates on the surface, and analyzes the changing pattern of this magnetic field to solve problems such as locating groundwater flow paths, detecting and locating underground pipelines, and exploring metal minerals.

[0003] The inversion technology of magnetoresistivity is an important method in geophysical exploration used to infer underground structure from surface observation data. Its core idea is to reconstruct the physical field distribution of underground resistivity, magnetic permeability, current density, wave velocity, etc. through the measurement results of known physical parameters combined with mathematical models and algorithms.

[0004] Progress in domestic magnetoresistivity inversion technology has been slow, limiting its engineering applications. Existing magnetoresistivity inversion techniques include: first implementing three-dimensional inversion of magnetoresistivity data using a regularized inversion approach, followed by the Gauss-Newton method to find the optimal solution for the objective function; implementing three-dimensional nonlinear inversion of borehole magnetoresistivity data using a nonlinear iterative least squares regularization technique and a Gauss-Newton algorithm to minimize the objective function; and implementing two nonlinear inversion techniques based on regularized inversion: the nonlinear conjugate gradient method and the limited-memory quasi-Newton method.

[0005] However, the inversion results of the magnetoresistivity method mentioned above have technical problems with poor accuracy, making it difficult to meet the requirements for fine imaging of underground conductors. The specific problems are as follows:

[0006] First, inversion is a typical underdetermined problem. Existing technologies typically use magnetic field data as the raw data for inversion imaging. This results in insufficient observational data, making it difficult to overcome the non-uniqueness of the inversion results. This leads to multiple solutions for the inversion imaging results (such as the location and morphology of dam leakage channels, metal pipelines, or ore bodies).

[0007] Second, existing technologies lack prior information constraints on inversion models, making it difficult to use known prior information to control inversion results, resulting in large deviations between the inversion imaging results (dam leakage channels, metal pipes or ore body locations and shapes) and the actual underground structure. Summary of the Invention

[0008] In response to the shortcomings of the existing technology, the present invention proposes a three-dimensional imaging method for underground conductors based on magnetoresistivity inversion technology. It not only improves the stability and reliability of the inversion imaging results by increasing the amount of data information and making full use of the information contained in multiple types of data; but also regulates the inversion imaging results by utilizing structural information such as the position and morphology of the field source and physical parameter information, thereby constraining the inversion of magnetic data to improve the stability and reliability of the calculation results, thereby solving the problem of three-dimensional fine imaging of underground conductors such as water and metal.

[0009] To achieve the above-mentioned object, the present invention designs a method for three-dimensional imaging of underground conductors based on magnetoresistivity inversion technology, which is particularly characterized by comprising the following steps:

[0010] S1) collecting on-site magnetic field three-directional component data at each observation point on the observation surface;

[0011] S2) preprocessing the three-directional component data of the magnetic field to extract the three-directional component data of the magnetic field for inversion;

[0012] S3) Establish the fitting difference function of the joint inversion data of the three-direction components of the magnetic field. The specific formula is as follows

[0013]

[0014] Where,

[0015] Φ d (m) represents the fitting difference function of the joint inversion data of the three directional components of the magnetic field,

[0016] m represents the current density to be solved in the inversion,

[0017] G i x (m) represents the forward kernel function of the x-direction component of the magnetic field,

[0018] G i y (m) represents the forward kernel function of the magnetic field y-direction component,

[0019] G i z (m) represents the forward kernel function of the z-direction component of the magnetic field,

[0020] d i obs-x Represents the collected magnetic field x-direction component data,

[0021] d i obs-y Represents the collected magnetic field y-direction component data,

[0022] d iobs-z Represents the collected magnetic field z-direction component data,

[0023] d max obs-k (k=x,y,z) represents the maximum value of the collected magnetic field x,y,z direction component data,

[0024] d min obs-k (k=x,y,z) represents the minimum value of the collected magnetic field x,y,z direction component data,

[0025] σ k (k=x,y,z) represents the data fitting error coefficient of the x,y,z direction components of the observation point. It is set to 3% to 5% by default according to the noise level of the observation data.

[0026] i represents the i-th surface magnetic field observation point,

[0027] N represents the number of magnetic field observation points;

[0028] S4) establishing a multi-source information constrained inversion model objective function, wherein the multi-source information constrained inversion model objective function constrains the inversion model by geological body spatial trend information, geological body spatial distribution information, geological body buried depth information, and geological body physical property information;

[0029] S5) constructing a Tikhonov regularized inversion objective function by combining the fitting difference function of the three-direction magnetic field component joint inversion data in step S3) with the multi-source information constraint inversion model objective function in step S4). The specific formula is as follows:

[0030] Φ(m)=Φ d (m)+βΦ m (m)

[0031] Where,

[0032] Φ(m) represents the Tikhonov regularized inversion objective function,

[0033] Φ d (m) represents the fitting difference function of the joint inversion data of the three directional components of the magnetic field,

[0034] Φ m (m) represents the objective function of the multi-source information constrained inversion model,

[0035] β is the regularization parameter, which indicates the relative weight of the two;

[0036] S6) solving the current density distribution corresponding to the minimum value of the Tikhonov regularized inversion objective function, and drawing a current density map of the geological body space. The area where the current density is higher than the overall value of the region is the spatial location of the underground conductor.

[0037] Furthermore, in S2), the preprocessing of the three-directional component data of the magnetic field includes interpolation, denoising, and field source separation.

[0038] Furthermore, in S4), the spatial trend constraint coefficient in the spatial trend information of the geological body is used to control and constrain the physical property differences allowed in a specific direction of the inversion model space as a whole;

[0039] Through the spatial structure constraint function in the geological body burial depth information, the physical property differences allowed in a specific area of the inversion model space are controlled and constrained locally;

[0040] By controlling the depth compensation function of the inverted current density depth distribution, the attenuation effect of the magnetic field with depth is overcome;

[0041] The physical property results of geological bodies obtained by geological analysis, drilling, and other geophysical exploration methods are converted into a physical property reference model for inversion, guiding the inversion in the direction determined by prior information.

[0042] Furthermore, in S4), the specific formula of the objective function of the multi-source information constrained inversion model is as follows:

[0043]

[0044] Where,

[0045] Φ m (m) represents the objective function of the multi-source information constrained inversion model,

[0046] m represents the current density to be solved in the inversion,

[0047] m ref Represents the physical property reference model used in the current inversion problem;

[0048] α k (k=s,x,y,z) represents the spatial direction constraint coefficient, where α s Corresponding to the minimum model objective function, α x ,α y ,α z Corresponding to the smooth model objective functions in the x, y, and z directions, respectively,

[0049] ω k (k=s,x,y,z) represents the spatial structure constraint function,

[0050] 1 / (z0+z)θ represents the depth compensation function, where z0 is the height of the observation point, z is the depth of the center of each grid unit, and θ is the attenuation coefficient related to the geometric shape of the underground conductor.

[0051] Furthermore, in S6), the minimum value of the Tikhonov regularized inversion objective function is iteratively solved using a nonlinear optimization method. In each iteration, a model current density change Δm is obtained by solving the problem, and the final current density result is m. n＝1 =m n +Δm, judge whether the final result of the current current density meets the inversion termination condition. If so, terminate; if not, set n=n+1 and return to the previous step.

[0052] Furthermore, in S6), Δm is solved by the following formula:

[0053]

[0054] Where,

[0055] J is the Jacobian matrix obtained by taking the first-order partial derivatives of the predicted magnetic field component data,

[0056] W d is a diagonal matrix, each element of which is the coefficient of the fitting difference function of the joint inversion data of the three directional components of the magnetic field in step S3),

[0057] C k (k=s,x,y,z)) is the spatial constraint weight coefficient matrix of the model space unit block (the jth one),

[0058] β is the regularization parameter,

[0059] d obs-k Represents the collected three-directional component data of the magnetic field.

[0060] Furthermore, in S6), the Jacobian matrix J and the diagonal matrix W d , spatial constraint weight coefficient matrix C k The specific formula is as follows:

[0061]

[0062] c k =α k *w k / (z0+z) θ

[0063] Where,

[0064] i represents the i-th surface magnetic field observation point,

[0065] j represents the jth underground inversion calculation grid,

[0066] d i pre-k (k=x,y,z) represents the predicted x,y,z direction component data,

[0067] m represents the current density to be solved in the inversion,

[0068] d max obs-k (k=x,y,z) represents the maximum value of the collected magnetic field x,y,z direction component data,

[0069] d min obs-k (k=x,y,z) represents the minimum value of the collected magnetic field x,y,z direction component data,

[0070] σ k (k=x,y,z) represents the data fitting error coefficient of the x,y,z direction components of the observation point. It is set to 3% to 5% by default according to the noise level of the observation data.

[0071] α k (k=s,x,y,z) represents the spatial direction constraint coefficient, where α s Corresponding to the minimum model objective function, α x ,α y ,α z Corresponding to the smooth model objective functions in the x, y, and z directions, respectively,

[0072] ω k (k=s,x,y,z) represents the spatial structure constraint function,

[0073] 1 / (z0+z) θ represents the depth compensation function, where z0 is the height of the observation point, z is the depth of the center of each grid unit, and θ is the attenuation coefficient related to the geometric shape of the underground conductor.

[0074] The advantages of the present invention are:

[0075] 1. Compared with the inversion calculation of single component data in the prior art, the present invention uses the three component data of x, y, and z of the observed magnetic field for inversion calculation at the same time, increasing the amount of data from N to 3N, thereby increasing the amount of data information and making full use of the information contained in the multi-type data to improve the stability and reliability of the inversion result, thus solving the problem of insufficient inversion data volume and insufficient solution reliability.

[0076] 2. The present invention uses the three-directional components of the magnetic field to construct a fitting difference function for the joint inversion data of the three-directional components of the magnetic field, establishes a joint inversion method for the three-directional components of the magnetic field data, increases the amount of data for magnetoresistivity inversion, and reduces the multi-solution problem of the inversion;

[0077] 3. Compared with the unconstrained inversion calculation in the prior art, the present invention proposes a constrained inversion method for three-dimensional magnetoresistivity, constructs a multi-source information constrained inversion model objective function, and uses structural information such as field source position and morphology and physical parameter information to regulate the inversion results, thereby constraining the magnetic data inversion to improve the stability and reliability of the calculation results. This overcomes the problem of a large amount of known information not being fully utilized, realizes the prior information constraint on the inversion process, and improves the inversion accuracy by regulating the inversion process.

[0078] The three-dimensional imaging method of underground conductors based on the magnetoresistivity inversion technology of the present invention improves the stability and reliability of the inversion imaging results, can be applied to the field of leakage detection in water conservancy projects, and solves the problem of three-dimensional fine positioning of leakage paths within the projects; at the same time, it can provide technical support for the monitoring of dam infiltration lines, effectively ensuring the safety of national water network project construction and operation; in addition, it can also be promoted and applied to groundwater structures, groundwater flow channels, underground pollution plumes, positioning of connected cracks or porous areas, monitoring of groundwater flow changes, monitoring of groundwater ion concentration changes, ground leaching solution monitoring and other groundwater and related structure detection, and has broad application prospects in many application fields such as urban underground engineering detection, environmental protection, disaster prevention and mitigation, etc. BRIEF DESCRIPTION OF THE DRAWINGS

[0079] Figure 1 is a flow chart of the present invention;

[0080] Figure 2 This is the model in Example 1 of the present invention;

[0081] Figure 3 The mesh division in the model of Example 1 of the present invention;

[0082] Figure 4 In order to use the method of the present invention Figure 2 Schematic diagram of the magnetic field in three directions for ground observations in the central model;

[0083] Figure 5 In order to use the method of the present invention Figure 2 Inversion imaging results of the inversion model;

[0084] Figure 6 It is the constraint effect of the spatial trend constraint coefficient on the inversion imaging results in the present invention;

[0085] Figure 7 The constraint effect of the spatial structure constraint function on the inversion imaging results in the present invention;

[0086] Figure 8 The constraint effect of the depth compensation function on the inversion imaging results in the present invention;

[0087] Figure 9 This is the constraint effect of the reference model on the inversion imaging results in the present invention. DETAILED DESCRIPTION

[0088] The present invention is further described in detail below with reference to the accompanying drawings and specific embodiments.

[0089] Example 1: Assume that there is a four-block water-bearing area underground. Figure 2 As shown, the water-bearing area of the four blocks is placed in a cubic low-resistance body of size 80*80*80m. The top surface of the cubic low-resistance body is buried at a depth of 80m and the conductivity is 0.1Sm -1 , the supply current is 1A, and the supply positions are (-600,0,0) and (600,0,0).

[0090] like Figure 1 As shown, the present invention provides a three-dimensional imaging method for underground conductors based on magnetoresistivity inversion technology, comprising the following steps:

[0091] S1) Collecting on-site the magnetic field three-directional component data Bx, By, and Bz of each observation point on the observation surface.

[0092] In this embodiment 1, Figure 3 As shown, there are three directions of magnetic field observed on the surface.

[0093] S2) Preprocessing the three-directional component data of the magnetic field to extract the three-directional component data of the magnetic field for inversion.

[0094] Specifically, the preprocessing of the three-directional component data of the magnetic field includes interpolation, denoising, and field source separation.

[0095] S3) Establish the fitting difference function of the joint inversion data of the three-direction components of the magnetic field. The specific formula is as follows

[0096]

[0097] Where,

[0098] Φ d (m) represents the fitting difference function of the joint inversion data of the three directional components of the magnetic field,

[0099] m represents the current density to be solved in the inversion,

[0100] G i x (m) represents the forward kernel function of the x-direction component of the magnetic field,

[0101] G i y (m) represents the forward kernel function of the magnetic field y-direction component,

[0102] G i z (m) represents the forward kernel function of the z-direction component of the magnetic field,

[0103] d i obs-x Represents the collected magnetic field x-direction component data,

[0104] d i obs-y Represents the collected magnetic field y-direction component data,

[0105] d i obs-z Represents the collected magnetic field z-direction component data,

[0106] d max obs-k (k=x,y,z) represents the maximum value of the collected magnetic field x,y,z direction component data,

[0107] d min obs-k (k=x,y,z) represents the minimum value of the collected magnetic field x,y,z direction component data,

[0108] σ k (k=x,y,z) represents the data fitting error coefficient of the x,y,z direction components of the observation point. It is set to 3% to 5% by default according to the noise level of the observation data.

[0109] i represents the i-th surface magnetic field observation point,

[0110] N represents the number of magnetic field observation points;

[0111] S4) establishing a multi-source information constrained inversion model objective function, wherein the multi-source information constrained inversion model objective function constrains the inversion model through geological body spatial trend information, geological body spatial distribution information, geological body burial depth information, and geological body physical property information.

[0112] Specifically, the spatial strike constraint coefficient in the spatial strike information of the geological body is used to control and constrain the physical property differences allowed in a specific direction of the inversion model space as a whole; the spatial structure constraint function in the burial depth information of the geological body is used to control and constrain the physical property differences allowed in a specific area of the inversion model space locally; the depth compensation function that controls the depth distribution of the inverted current density is used to overcome the attenuation effect of the magnetic field with depth; the physical property results of the geological body obtained by geological analysis, drilling, and other geophysical detection methods are converted into a physical property reference model for inversion, guiding the inversion in the direction determined by the prior information.

[0113] Specifically, the specific formula of the objective function of the multi-source information constraint inversion model is as follows:

[0114]

[0115] Where,

[0116] Φ m (m) represents the objective function of the multi-source information constrained inversion model,

[0117] m represents the current density to be solved in the inversion,

[0118] m ref Represents the physical property reference model used in the current inversion problem;

[0119] α k (k=s,x,y,z) represents the spatial direction constraint coefficient, where α s Corresponding to the minimum model objective function, α x ,α y ,α z Corresponding to the smooth model objective functions in the x, y, and z directions, respectively,

[0120] ω k (k=s,x,y,z) represents the spatial structure constraint function,

[0121] 1 / (z0+z) θ represents the depth compensation function, where z0 is the height of the observation point, z is the depth of the center of each grid unit, and θ is the attenuation coefficient related to the geometric shape of the underground conductor.

[0122] Specifically, the spatial orientation constraint coefficient α k The weight coefficients of each item in the control model objective function are used to control the physical property differences allowed in a specific direction of the model space. Specifically, through α k The value of the (k=s,x,y,z) coefficient is used to balance the relative weight between the minimum model objective function and the smoothest model objective function. s ,α x ,αy ,α z The default value is 1,1,1,1, and the range is 0.001 to 1000. The smaller the value, the greater the weight. When the spatial orientation of the geological body is known based on external data, this directional coefficient can be reduced to increase the weight of the allowable physical property changes in that direction, thus incorporating the spatial orientation constraint of the geological body into the inversion.

[0123] Specifically, the spatial structure constraint function ω k (k=s,x,y,z) is a flexible means of inversion spatial structure constraint, which can locally control the physical property differences allowed in a specific area of the model space. It can be set independently for each grid cell, aiming to adjust the degree of difference between each local grid cell and the physical property reference model, thereby improving the accuracy of the inversion results. Specifically, by converting the geometric parameter information describing the field source body, such as the known boundary position, top surface depth, and inclination, into a spatial weighting function, ω s ,ω x ,ω y ,ω z The default value is 1,1,1,1, and the range is 0.001 to 1000. The smaller the value, the greater the weight. When the spatial distribution area of the geological body is known based on external data, the regional coefficient can be reduced to increase the weight of the allowable physical property changes in the area, and the spatial geometric structure constraints of the geological body can be added to the inversion.

[0124] Specifically, θ in the depth compensation function is an attenuation coefficient related to the geometric shape of the underground conductor, and a value of 3 is recommended.

[0125] Specifically, the physical property reference model m ref You can set a specific physical property value m0 for each unit independently, or you can limit the range of physical property changes of each grid m min <m0<m max By converting the physical property results obtained by geological analysis, drilling, and other geophysical exploration methods into a physical property reference model for inversion, the inversion is guided in the direction determined by the prior information, which plays a role in suppressing the multi-solution problem of inversion.

[0126] S5) constructing a Tikhonov regularized inversion objective function by combining the fitting difference function of the three-direction magnetic field component joint inversion data in step S3) with the multi-source information constraint inversion model objective function in step S4). The specific formula is as follows:

[0127] Φ(m)=Φ d (m)+βΦ m (m)

[0128] Where,

[0129] Φ(m) represents the Tikhonov regularized inversion objective function,

[0130] Φ d (m) represents the fitting difference function of the joint inversion data of the three directional components of the magnetic field,

[0131] Φ m (m) represents the objective function of the multi-source information constrained inversion model,

[0132] β is the regularization parameter, which indicates the relative weight of the two.

[0133] Specifically, in the Tikhonov regularized inversion objective function, the former establishes the relationship between the observed data and the physical property parameters in the model space through the forward formula, and the latter serves as a regularized stability term and can control the distribution characteristics of the physical properties in the model space; the relative weight of the two is represented by the regularization parameter β.

[0134] S6) Solving the current density distribution corresponding to the minimum value of the Tikhonov regularized inversion objective function and plotting a current density map of the geological volume. Regions where the current density is higher than the overall regional value are the spatial locations of underground conductors, specifically, underground dam leakage channels, metal pipelines, or ore bodies.

[0135] Specifically, the minimum value of the Tikhonov regularized inversion objective function is iteratively solved using a nonlinear optimization method. In each iteration, a model current density change Δm is obtained by solving the problem, and then the final current density result is m. n＝1 =m n +Δm, judge whether the final result of the current current density meets the inversion termination condition. If so, terminate; if not, set n=n+1 and return to the previous step.

[0136] Specifically, Δm is solved by the following formula:

[0137]

[0138] Where,

[0139] J is the Jacobian matrix obtained by taking the first-order partial derivatives of the predicted magnetic field component data,

[0140] W d is a diagonal matrix, each element of which is the coefficient of the fitting difference function of the joint inversion data of the three directional components of the magnetic field in step S3),

[0141] C k (k=s,x,y,z)) is the spatial constraint weight coefficient matrix of the model space unit block (the jth one),

[0142] β is the regularization parameter,

[0143] d obs-k Represents the collected three-directional component data of the magnetic field.

[0144] Specifically, the Jacobian matrix J, the diagonal matrix W d , spatial constraint weight coefficient matrix C k The specific formula is as follows:

[0145]

[0146] c k =α k *w k / (z0+z) θ

[0147] Where,

[0148] i represents the i-th surface magnetic field observation point,

[0149] j represents the jth underground inversion calculation grid,

[0150] d i pre-k (k=x,y,z) represents the predicted x,y,z direction component data,

[0151] m represents the current density to be solved in the inversion,

[0152] d max obs-k (k=x,y,z) represents the maximum value of the collected magnetic field x,y,z direction component data,

[0153] d min obs-k (k=x,y,z) represents the minimum value of the collected magnetic field x,y,z direction component data,

[0154] α k (k=s,x,y,z) represents the spatial direction constraint coefficient, where α s Corresponding to the minimum model objective function, α x ,α y ,α z Corresponding to the smooth model objective functions in the x, y, and z directions, respectively,

[0155] σ k (k=x,y,z) represents the data fitting error coefficient of the x,y,z direction components of the observation point. It is set to 3% to 5% by default according to the noise level of the observation data.

[0156] ω k (k=s,x,y,z) represents the spatial structure constraint function,

[0157] 1 / (z0+z) θ represents the depth compensation function, where z0 is the height of the observation point, z is the depth of the center of each grid unit, and θ is the attenuation coefficient related to the geometric shape of the underground conductor.

[0158] In this embodiment 1, the inversion imaging result obtained by using the present invention is as follows: Figure 4 shown. Figure 4 This is a longitudinal section of a cubic low-resistance body. The underground current density distribution is calculated through inversion. The red area in the figure is the area with the maximum current density, which is the speculated position of the conductor, thus realizing the positioning of the leakage area inside the underground pipeline.

[0159] The following discusses the constraint effects of the spatial orientation constraint coefficient, spatial structure constraint function, depth compensation function, and physical property reference model on the inversion results in the present invention.

[0160] Example 2

[0161] like Figure 6 As shown, it is the spatial direction constraint coefficient α in the present invention k The constraint effect of (k=s,x,y,z) on the inversion results.

[0162] This embodiment sets up a combined model of three cubes, whose spatial positions are shown in black squares, keeping α s ,α x ,α y The default value is 1, set α z Different values of 1, 10, and 100 gradually control the inversion calculated current to become narrower in the horizontal direction and wider in the vertical direction. This embodiment shows that the inversion calculation results can be controlled to be distributed in a specified direction through the spatial direction constraint coefficient.

[0163] Example 3

[0164] like Figure 7 As shown, it is the spatial structure constraint function ω in the present invention. k The constraint effect of (k=s,x,y,z) on the inversion results.

[0165] This example assumes the existence of a known geological interface. The high conductivity areas interpreted from the drillhole data are primarily distributed below this interface, as shown by the black dashed line in the figure. Using this interface as the dividing line, different spatial structure constraint coefficients are set. The high current density in the inversion results is distributed below this interface, enabling regional control of the inversion results. The calculated results meet known prior geological information and are more reliable.

[0166] Example 4

[0167] like Figure 8As shown, it is the depth compensation function 1 / (z0+z) in the present invention. θ Constraint effects on inversion results.

[0168] The depth compensation coefficient θ = 1, 2, 3, or 4 varies. As the depth compensation coefficient increases, the high current density area calculated by the inversion calculation gradually shifts downward. θ = 1 or 2 provide insufficient depth compensation, while θ = 4 overcompensates. When θ = 3, the calculated high current density area best corresponds to the theoretical model depth position. This example demonstrates that the depth compensation function can be used to control the depth distribution of the inversion calculation results, and the empirically optimal value θ = 3 is obtained.

[0169] Example 5

[0170] like Figure 9 As shown, this is the physical property reference model m in the present invention. ref Constraint effects on inversion results.

[0171] This example assumes that there are already results obtained by other geophysical methods. The results obtained by other geophysical methods are converted into a physical property reference model. After the physical property reference model is added to the inversion calculation, the inverted current density distribution in the blue low current density area is similar to the reference model. The physical property reference model has a good constraint on the inversion results. In addition, there is a significant difference between the red high current density area and the reference model, and the inversion results are not completely dependent on the reference model. Physical property reference model m ref It effectively constrains the inversion results, but does not force the inversion results to be similar.

[0172] The three-dimensional imaging method of underground conductors based on the magnetoresistivity inversion technology of the present invention improves the stability and reliability of the inversion imaging results, and is applied to the field of leakage detection in water conservancy projects to solve the problem of three-dimensional precise positioning of leakage paths within the projects; at the same time, it can provide technical support for the monitoring of dam infiltration lines, effectively ensuring the safety of national water network project construction and operation; in addition, it can also be promoted and applied to groundwater structures, groundwater flow channels, underground pollution plumes, positioning of connected cracks or porous areas, monitoring of groundwater flow changes, monitoring of groundwater ion concentration changes, ground leaching solution monitoring and other groundwater and related structure detection, and has broad application prospects in many application fields such as urban underground engineering detection, environmental protection, disaster prevention and mitigation, etc.

[0173] The above embodiments are preferred implementation modes of the present invention, but the implementation modes of the present invention are not limited to the above embodiments. Any other changes, modifications, substitutions, combinations, and simplifications that do not deviate from the spirit and principles of the present invention should be considered as equivalent replacement methods and are included in the scope of protection of the present invention.

Claims

1. A method for three-dimensional imaging of underground conductors based on magnetoresistivity inversion technology, characterized in that: The steps include: S1) collecting on-site magnetic field three-directional component data at each observation point on the observation surface; S2) preprocessing the three-directional component data of the magnetic field to extract the three-directional component data of the magnetic field for inversion; S3) Establish the fitting difference function of the joint inversion data of the three-direction components of the magnetic field. The specific formula is as follows Where, Φ d (m) represents the fitting difference function of the joint inversion data of the three directional components of the magnetic field, m represents the current density to be solved in the inversion, G i x (m) represents the forward kernel function of the x-direction component of the magnetic field, G i y (m) represents the forward kernel function of the magnetic field y-direction component, G i z (m) represents the forward kernel function of the z-direction component of the magnetic field, d i obs-x Represents the collected magnetic field x-direction component data, d i obs-y Represents the collected magnetic field y-direction component data, d i obs-z Represents the collected magnetic field z-direction component data, d max obs-k (k=x,y,z) represents the maximum value of the collected magnetic field x,y,z direction component data, d min obs-k (k=x,y,z) represents the minimum value of the collected magnetic field x,y,z direction component data, σ k (k=x,y,z) represents the data fitting error coefficient of the x,y,z direction components of the observation point. It is set to 3% to 5% by default according to the noise level of the observation data. i represents the i-th surface magnetic field observation point, N represents the number of magnetic field observation points; S4) establishing a multi-source information constrained inversion model objective function, wherein the multi-source information constrained inversion model objective function constrains the inversion model by geological body spatial trend information, geological body spatial distribution information, geological body buried depth information, and geological body physical property information; S5) constructing a Tikhonov regularized inversion objective function by combining the fitting difference function of the three-direction magnetic field component joint inversion data in step S3) with the multi-source information constraint inversion model objective function in step S4). The specific formula is as follows: Φ(m)=Φ d (m)+βΦ m (m) Where, Φ(m) represents the Tikhonov regularized inversion objective function, Φ d (m) represents the fitting difference function of the joint inversion data of the three directional components of the magnetic field, Φ m (m) represents the objective function of the multi-source information constrained inversion model, β is the regularization parameter, which indicates the relative weight of the two; S6) solving the current density distribution corresponding to the minimum value of the Tikhonov regularized inversion objective function, and drawing a current density map of the geological body space. The area where the current density is higher than the overall value of the region is the spatial location of the underground conductor.

2. The method for three-dimensional imaging of underground conductors based on magnetoresistivity inversion technology according to claim 1, characterized in that: In S2), the preprocessing of the three-directional component data of the magnetic field includes interpolation, denoising, and field source separation.

3. The method for three-dimensional imaging of underground conductors based on magnetoresistivity inversion technology according to claim 1, characterized in that: In S4), the spatial trend constraint coefficient in the spatial trend information of the geological body is used to control and constrain the physical property differences allowed in a specific direction of the inversion model space as a whole; Through the spatial structure constraint function in the geological body burial depth information, the physical property differences allowed in a specific area of the inversion model space are controlled and constrained locally; By controlling the depth compensation function of the inverted current density depth distribution, the attenuation effect of the magnetic field with depth is overcome; The physical property results of geological bodies obtained by geological analysis, drilling, and other geophysical exploration methods are converted into a physical property reference model for inversion, guiding the inversion in the direction determined by prior information.

4. The method for three-dimensional imaging of underground conductors based on magnetoresistivity inversion technology according to claim 3, characterized in that: In S4), the specific formula of the objective function of the multi-source information constraint inversion model is as follows: Where, Φ m (m) represents the objective function of the multi-source information constrained inversion model, m represents the current density to be solved in the inversion, m ref Represents the physical property reference model used in the current inversion problem; α k (k=s,x,y,z) represents the spatial direction constraint coefficient, where α s Corresponding to the minimum model objective function, α x ,α y ,α z Corresponding to the smooth model objective functions in the x, y, and z directions, respectively, ω k (k=s,x,y,z) represents the spatial structure constraint function, 1 / (z0+z) θ represents the depth compensation function, where z0 is the height of the observation point, z is the depth of the center of each grid unit, and θ is the attenuation coefficient related to the geometric shape of the underground conductor.

5. The method for three-dimensional imaging of underground conductors based on magnetoresistivity inversion technology according to claim 1, characterized in that: In S6), the minimum current density of the Tikhonov regularized inversion objective function is iteratively solved using a nonlinear optimization method. In each iteration, a model current density change Δm is obtained by solving the problem, and the final current density result is m. n＝1 =m n +Δm, judge whether the final result of the current current density meets the inversion termination condition. If so, terminate; if not, set n=n+1 and return to the previous step.

6. The method for three-dimensional imaging of underground conductors based on magnetoresistivity inversion technology according to claim 5, characterized in that: In S6), Δm is solved by the following formula: Where, J is the Jacobian matrix obtained by taking the first-order partial derivatives of the three-directional component data of the predicted magnetic field, W d is a diagonal matrix, each element of which is the coefficient of the fitting difference function of the joint inversion data of the three directional components of the magnetic field in step S3), C k (k=s,x,y,z)) is the spatial constraint weight coefficient matrix of the model space unit block (the jth one), β is the regularization parameter, d obs-k Represents the collected three-directional component data of the magnetic field.

7. The method for three-dimensional imaging of underground conductors based on magnetoresistivity inversion technology according to claim 6, characterized in that: In S6), the Jacobian matrix J and the diagonal matrix W d , spatial constraint weight coefficient matrix C k The specific formula is as follows: c k =α k *In k / (z0+z) θ Where, i represents the i-th surface magnetic field observation point, j represents the jth underground inversion calculation grid, d i pre-k (k=x,y,z) represents the predicted x,y,z direction component data, m represents the current density to be solved in the inversion, d max obs-k (k=x,y,z) represents the maximum value of the collected magnetic field x,y,z direction component data, d min obs-k (k=x,y,z) represents the minimum value of the collected magnetic field x,y,z direction component data, σ k (k=x,y,z) represents the data fitting error coefficient of the x,y,z direction components of the observation point. It is set to 3% to 5% by default according to the noise level of the observation data. α k (k=s,x,y,z) represents the spatial direction constraint coefficient, where α s Corresponding to the minimum model objective function, α x ,α y ,α z Corresponding to the smooth model objective functions in the x, y, and z directions, respectively, ω k (k=s,x,y,z) represents the spatial structure constraint function, 1 / (z0+z) θ represents the depth compensation function, where z0 is the height of the observation point, z is the depth of the center of each grid unit, and θ is the attenuation coefficient related to the geometric shape of the underground conductor.

Images