A joint inversion method and related equipment for alternating iterative inversion of cross-aperture radar travel time CT and high-density resistivity.

Abstract

This invention discloses a joint inversion method and related equipment for alternating iterative inversion using cross-aperture radar travel-time CT and high-density resistivity methods: acquiring electromagnetic wave data and apparent resistivity data within the target observation area; constructing the objective function for the joint inversion using alternating iterative inversion of cross-aperture radar travel-time CT and high-density resistivity methods; establishing initial velocity and initial resistivity models; acquiring target forward response data and the target Jacobian matrix; acquiring the cross-gradient partial derivatives of the overlapping observation area; establishing a large matrix and acquiring the update amounts of the velocity and resistivity models; acquiring the current velocity and current resistivity models; returning to the steps of acquiring target forward response data using cross-aperture radar travel-time CT and high-density resistivity methods until the number of iterations reaches its maximum value or the fitting difference between the target forward response data and the observation data reaches a preset threshold, and outputting the inversion result. This invention can improve the accuracy of alternating iterative joint inversion and can be widely applied in the field of geophysical exploration technology.

Description

Technical Field

[0001] This invention relates to the field of geophysical exploration technology, and in particular to a method and related equipment for joint inversion of cross-hole radar travel-time CT and high-density resistivity alternating iterative methods. Background Technology

[0002] Geophysical exploration is a crucial method for probing the structure, physical properties, and location of subsurface media using geophysical methods. Inversion is a core step in geophysical data processing, aiming to quantitatively reconstruct the physical parameters of subsurface media from observational data, providing a quantitative basis for subsequent geological interpretation. Since the ambiguity of inversion results is one of the most fundamental challenges in geophysical exploration, and single inversion methods cannot effectively reduce this ambiguity, cross-gradient joint inversion is an important means of integrating geophysical methods for interpretation. However, synchronous iterative cross-gradient joint inversion suffers from problems such as excessively large joint inversion matrix size, unstable solution process, and difficulty in optimizing method weighting factors. Step-by-step iterative cross-gradient joint inversion has drawbacks, including difficulty in balancing the weights of data fitting terms and structural constraint terms, and it cannot effectively reduce the ambiguity of inversion results. Summary of the Invention

[0003] In view of this, the main objective of the embodiments of the present invention is to provide a method and related equipment for joint inversion of cross-aperture radar travel time CT and high-density resistivity alternating iterative method, in order to solve at least one of the problems of the prior art. The present invention can improve the accuracy of joint inversion of alternating iterative cross-gradient method.

[0004] To achieve the above objectives, one aspect of the present invention provides a joint inversion method for alternating iterative inversion of transap-scan radar travel time CT and high-density resistivity, the method comprising: Acquire electromagnetic wave data and apparent resistivity data within the target observation area; Based on the electromagnetic wave data and the apparent resistivity data, an objective function is constructed by alternating iterative joint inversion of cross-aperture radar travel-time CT and high-density resistivity method. Establish an initial velocity model and an initial resistivity model; use the initial velocity model and the initial resistivity model as models to be processed; Based on the model to be processed, obtain the target forward modeling response data of the cross-aperture radar travel time CT and the high-density resistivity method, and based on the target forward modeling response data, obtain the target Jacobian matrix of the cross-aperture radar travel time CT and the high-density resistivity method. By minimizing the objective function, the cross gradient partial derivatives of the overlapping observation regions within the target observation region are obtained; Based on the target forward response data, the target Jacobian matrix, and the cross gradient partial derivatives, a large matrix is ​​established for the cross-aperture radar travel time CT and the high-density resistivity method, and the velocity model update and resistivity model update are obtained based on the large matrix. Based on the model to be processed, the velocity model update amount, and the resistivity model update amount, obtain the current velocity model and the current resistivity model; The current velocity model and the current resistivity model are used as the model to be processed. The process of obtaining the target forward modeling response data of the cross-aperture radar travel time CT and the high-density resistivity method based on the model to be processed is repeated until the number of iterations reaches the maximum value or the fitting difference between the target forward modeling response data and the observation data reaches a preset threshold. The joint inversion result of alternating iterations is then output.

[0005] In some embodiments, acquiring electromagnetic wave data and apparent resistivity data within the target observation area includes the following steps: The apparent resistivity data is obtained based on the high-density resistivity survey lines pre-arranged in the target observation area of ​​the underground karst. Two boreholes are drilled in the target observation area, and a cross-hole radar receiving antenna and a cross-hole radar transmitting antenna are respectively installed in the two boreholes to form the overlapping observation area; Based on the cross-hole radar receiving antenna and the cross-hole radar transmitting antenna in the overlapping observation area of ​​the underground karst, the electromagnetic wave data is collected by moving along the depth direction of the borehole at a fixed measuring point spacing. The distance between the measuring points is less than half the wavelength of the main frequency of the cross-hole radar.

[0006] In some embodiments, the step of constructing an objective function jointly inverted by alternating iterative methods of transaperture radar travel-time CT and high-density resistivity method based on the electromagnetic wave data and the apparent resistivity data includes the following steps: Based on the model gradient within the overlapping observation area using cross-aperture radar travel-time CT and high-density resistivity method, local structural constraints are constructed. Acquire the observation data from the electromagnetic wave data and the apparent resistivity data, including the cross-aperture radar travel-time CT and the high-density resistivity method. Obtain the reciprocal matrix of the standard deviation of the observation data to get the data weighting matrix of cross-aperture radar travel time CT and high-density resistivity method; The discretized inversion grid structure of the target observation area is obtained. Based on the discretized inversion grid structure, a second-order difference operator is constructed to obtain the model weighting matrix of the cross-aperture radar travel time CT and the high-density resistivity method. Based on the local structural constraints, the observation data, the data weighting matrix, and the model weighting matrix, the objective function is constructed by alternating iterations of the cross-aperture radar travel-time CT and the high-density resistivity method under local structural constraints.

[0007] In some embodiments, establishing the initial velocity model and the initial resistivity model includes the following steps: Obtain the average velocity within the overlapping observation area, and set the average velocity as the initial velocity model; Obtain the average resistivity within the target observation area, and set the average resistivity as the initial resistivity model.

[0008] In some embodiments, obtaining target forward modeling response data based on the cross-aperture radar travel-time CT and high-density resistivity method according to the model to be processed includes the following steps: The finite difference method is used to solve the equation of the process function, and the forward modeling response of the cross-aperture radar travel-time CT is calculated to obtain the first forward modeling response data of the cross-aperture radar travel-time CT. The finite element method is used to perform forward response calculations on the high-density resistivity method to obtain the second forward response data of the high-density resistivity method; The target forward modeling response data includes the first forward modeling response data and the second forward modeling response data.

[0009] In some embodiments, obtaining the target Jacobian matrix using the cross-aperture radar travel-time CT and the high-density resistivity method based on the target forward response data includes the following steps: The first Jacobian matrix of the cross-aperture radar travel time CT is obtained by using the interchange theorem; The second Jacobian matrix of the high-density resistivity method is obtained by using the interchange theorem. The target Jacobian matrix includes the first Jacobian matrix and the second Jacobian matrix.

[0010] In some embodiments, establishing a large matrix for cross-aperture radar travel-time CT and high-density resistivity method based on the target forward response data, the target Jacobian matrix, and the cross gradient partial derivatives includes the following steps: Acquire the observation data from the electromagnetic wave data and the apparent resistivity data, including the cross-aperture radar travel-time CT and the high-density resistivity method. Obtain the reciprocal matrix of the standard deviation of the observation data to get the data weighting matrix of cross-aperture radar travel time CT and high-density resistivity method; Obtain the residual between the observed data and the target forward response data; The discretized inversion grid structure of the target observation area is obtained. Based on the discretized inversion grid structure, a second-order difference operator is constructed to obtain the model weighting matrix of the cross-aperture radar travel time CT and the high-density resistivity method. The large matrix is ​​constructed based on the data weighting matrix, the residual, the model weighting matrix, the target Jacobian matrix, and the cross gradient partial derivatives.

[0011] To achieve the above objectives, another aspect of the present invention proposes a joint inversion device for alternating iterative inversion of cross-aperture radar travel time CT and high-density resistivity, the device comprising: The data acquisition module is used to acquire electromagnetic wave data and apparent resistivity data within the target observation area; The objective function construction module is used to construct an objective function jointly inverted by alternating iterations of cross-aperture radar travel-time CT and high-density resistivity method based on the electromagnetic wave data and the apparent resistivity data. An initial model building module is used to establish an initial velocity model and an initial resistivity model; the initial velocity model and the initial resistivity model are used as models to be processed. The forward modeling module is used to obtain the target forward modeling response data of the cross-aperture radar travel time CT and the high-density resistivity method based on the model to be processed, and to obtain the target Jacobian matrix of the cross-aperture radar travel time CT and the high-density resistivity method based on the target forward modeling response data. The cross-gradient partial derivative acquisition module is used to obtain the cross-gradient partial derivatives of the overlapping observation regions within the target observation region by minimizing the objective function; The model update quantity acquisition module is used to establish a large matrix of cross-aperture radar travel time CT and high-density resistivity method based on the target forward response data, the target Jacobian matrix and the cross gradient partial derivatives, and to acquire the velocity model update quantity and resistivity model update quantity based on the large matrix. The model update module is used to obtain the current velocity model and the current resistivity model based on the model to be processed, the velocity model update amount, and the resistivity model update amount. The iterative inversion module is used to take the current velocity model and the current resistivity model as the model to be processed, and return the step of obtaining the target forward modeling response data of the cross-aperture radar travel time CT and the high-density resistivity method according to the model to be processed, until the number of iterations reaches the maximum value, or the fitting difference between the target forward modeling response data and the observation data reaches a preset threshold, and output the joint inversion result of alternating iterations.

[0012] To achieve the above objectives, another aspect of the present invention provides an electronic device, the electronic device including a memory and a processor, the memory storing a computer program, and the processor executing the computer program to implement the method described above.

[0013] To achieve the above objectives, another aspect of the present invention provides a computer-readable storage medium storing a computer program that, when executed by a processor, implements the methods described above.

[0014] To achieve the above objectives, another aspect of the present invention provides a computer program product or computer program that includes computer instructions stored in a computer-readable storage medium. A processor of a computer device can read the computer instructions from the computer-readable storage medium and execute the computer instructions to cause the computer device to perform the aforementioned method.

[0015] The embodiments of the present invention include at least the following beneficial effects: The present invention provides a method and related equipment for joint inversion of cross-aperture radar travel-time CT and high-density resistivity alternating iteration. This scheme obtains electromagnetic wave data and apparent resistivity data within the target observation area, providing a multi-source physical field observation data foundation for joint inversion; based on electromagnetic wave data and apparent resistivity data, an objective function for joint inversion of cross-aperture radar travel-time CT and high-density resistivity method alternating iteration is constructed. By introducing an alternating iteration strategy, the ambiguity of the inversion results of a single method is effectively reduced without relying on complex rock physical property relationships; an initial velocity model and an initial resistivity model are established, and a uniform half-space initialization model is adopted to simplify the inversion starting conditions and provide a stable calculation benchmark for subsequent cross-gradient structure coupling; the initial velocity model and the initial resistivity model are used as models to be processed; based on the models to be processed, the target forward modeling response data of cross-aperture radar travel-time CT and high-density resistivity method are obtained, and the target Jacobian matrix of cross-aperture radar travel-time CT and high-density resistivity method is obtained based on the target forward modeling response data. The Jacobian matrix is ​​solved efficiently using the interchange theorem, significantly reducing large-scale matrix operations. This approach reduces computational complexity and improves inversion efficiency. It minimizes the objective function to obtain the cross-gradient partial derivatives of overlapping observation regions within the target observation area. Based on the target forward response data, the target Jacobian matrix, and the cross-gradient partial derivatives, it establishes large matrices for the cross-aperture radar travel-time CT and high-density resistivity method, and obtains the velocity model update and resistivity model update based on these large matrices. In the objective function, it sets an adaptive regularization factor for the structural constraint term to dynamically balance the weights of the data fitting term and the structural constraint term, and sets a cross-gradient weight factor for the cross-gradient term to dynamically balance the constraint effects of the velocity and resistivity iterative models, thus improving the reliability of the inversion results. Based on the model to be processed, the velocity model update, and the resistivity model update, it obtains the current velocity model and the current resistivity model. Using the current velocity model and the current resistivity model as the model to be processed, it returns to the steps of obtaining the target forward response data for the cross-aperture radar travel-time CT and high-density resistivity method based on the model to be processed, until the number of iterations reaches its maximum value or the fitting difference between the target forward response data and the observation data reaches a preset threshold, outputting the joint inversion results of alternating iterations. This invention effectively avoids erroneous inversion results caused by the direct coupling of physical property parameters with large differences in magnitude by normalizing the partial derivatives of the cross gradient, thereby improving the accuracy of the joint inversion of alternating iterative cross gradients. Attached Figure Description

[0016] To more clearly illustrate the technical solutions in the embodiments of the present invention, the accompanying drawings used in the description of the embodiments will be briefly introduced below. Obviously, the accompanying drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0017] Figure 1 This is a flowchart of the cross-aperture radar travel time CT and high-density resistivity alternating iterative joint inversion method provided in this embodiment of the invention; Figure 2 This is a schematic diagram of the cross gradient constraint within the overlapping observation area provided in an embodiment of the present invention; Figure 3 This is a theoretical model diagram of resistivity and velocity provided in an embodiment of the present invention; Figure 4 This is a flowchart illustrating the specific implementation steps of the cross-aperture radar travel time CT and high-density resistivity alternating iterative joint inversion provided in this embodiment of the invention; Figure 5 This is a schematic diagram of the individual inversion results provided in an embodiment of the present invention; Figure 6 This is a schematic diagram of the joint inversion results provided in an embodiment of the present invention; Figure 7 This is a schematic diagram of the hardware structure of the electronic device provided in an embodiment of the present invention. Detailed Implementation

[0018] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in 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. The embodiments described in the following exemplary embodiments do not represent all embodiments consistent with those of this invention; they are merely examples of apparatuses and methods consistent with some aspects of the embodiments of this invention as detailed in the appended claims.

[0019] It should be noted that although functional modules are divided in the system diagram and a logical order is shown in the flowchart, in some cases, the steps shown or described may be performed in a different order than the module division in the system or the order in the flowchart. The terms "first / S100" and "second / S200" in the specification, claims, and the foregoing drawings may be used herein to describe various concepts, but unless specifically stated otherwise, these concepts are not limited by these terms. These terms are used only to distinguish one concept from another. For example, first information may also be referred to as second information without departing from the scope of the embodiments of the invention, and similarly, second information may also be referred to as first information. Depending on the context, the words "if" or "when" as used herein may be interpreted as "when," "in response to a determination," or "in the event of a determination."

[0020] The terms “at least one,” “multiple,” “each,” “any,” etc., used in this invention, “at least one” includes one, two, or more than two; “multiple” includes two or more than two; “each” refers to each of the corresponding multiple; and “any” refers to any one of the multiple.

[0021] Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this invention pertains. The terminology used herein is for the purpose of describing embodiments of the invention only and is not intended to limit the scope of the invention.

[0022] Before providing a detailed description of the embodiments of the present invention, some of the nouns and terms involved in the embodiments of the present invention will be explained first. The nouns and terms involved in the embodiments of the present invention are subject to the following interpretations.

[0023] The high-density resistivity method is a geophysical exploration method that uses a multi-electrode system to measure the resistivity distribution of underground soil and rock at high resolution. It is used to quickly and accurately obtain electrical structure information of underground soil and rock masses.

[0024] Underground karst is an underground cave system formed by the chemical and physical action of water on soluble rocks (such as limestone), mainly including caves and underground rivers.

[0025] The Cheng-Jun equation is a class of nonlinear partial differential equations that describe the propagation of waves in a medium. It is used to connect wave optics and geometric optics and is widely used in seismic exploration and physical optics.

[0026] Among the related technologies, there are single inversion methods and cross-gradient joint inversion methods. Among them, the single inversion method cannot effectively reduce the ambiguity of the inversion results, while the existing synchronous iterative cross-gradient joint inversion has an excessively large joint inversion matrix, unstable solution process, and difficulty in optimizing the selection of method weight factors. The step-by-step iterative cross-gradient joint inversion has difficulty balancing the weights of data fitting terms and structural constraint terms, and cannot effectively reduce the ambiguity of the inversion results.

[0027] In view of this, this invention provides a method and related equipment for joint inversion of cross-aperture radar travel-time CT and high-density resistivity using alternating iterative methods. This approach introduces an alternating iterative joint inversion strategy into the joint inversion method of cross-aperture radar travel-time CT and high-density resistivity with local structural constraints. It does not rely on complex rock property relationships and can effectively reduce the ambiguity of the inversion results. Secondly, an adaptive regularization factor is set for the structural constraint term to balance the data fitting term and the structural constraint term; a cross-gradient weight factor is set for the cross-gradient term to balance the constraint effects of the velocity and resistivity iterative models, thereby improving the reliability of the inversion results. Simultaneously, normalizing the cross-gradient function and partial derivatives can effectively avoid erroneous inversion results caused by the direct coupling of physical property parameters with large differences in magnitude.

[0028] The cross-aperture radar travel-time CT and high-density resistivity alternating iterative joint inversion method provided in this invention relates to the field of geophysical exploration technology. This method can be applied to terminals, servers, or software running on either. In some embodiments, the terminal can be a smartphone, tablet, laptop, desktop computer, smart speaker, smartwatch, or vehicle terminal, but is not limited to these. The server can be configured as an independent physical server, a server cluster or distributed system composed of multiple physical servers, or a cloud server providing basic cloud computing services such as cloud services, cloud databases, cloud computing, cloud functions, cloud storage, network services, cloud communication, middleware services, domain name services, security services, CDN, and big data and artificial intelligence platforms. The server can also be a node server in a blockchain network. The software can be an application implementing the cross-aperture radar travel-time CT and high-density resistivity alternating iterative joint inversion method, but is not limited to the above forms.

[0029] Figure 1 This is an optional flowchart of a joint inversion method for alternating iterative inversion of cross-aperture radar travel time CT and high-density resistivity provided in an embodiment of the present invention. Figure 1 The method may include, but is not limited to, steps S100 to S800: Step S100: Obtain electromagnetic wave data and apparent resistivity data within the target observation area; Step S200: Based on electromagnetic wave data and apparent resistivity data, construct the objective function for joint inversion by alternating iterations of cross-aperture radar travel time CT and high-density resistivity method; Step S300: Establish the initial velocity model and the initial resistivity model; use the initial velocity model and the initial resistivity model as the models to be processed; Step S400: Based on the model to be processed, obtain the cross-aperture radar travel time CT and the target forward modeling response data of the high-density resistivity method, and based on the target forward modeling response data, obtain the cross-aperture radar travel time CT and the target Jacobian matrix of the high-density resistivity method. Step S500: By minimizing the objective function, obtain the cross gradient partial derivatives of the overlapping observation regions within the target observation region; Step S600: Based on the target forward response data, the target Jacobian matrix, and the cross gradient partial derivatives, establish a large matrix for the cross-aperture radar travel time CT and the high-density resistivity method, and obtain the velocity model update and resistivity model update based on the large matrix. Step S700: Based on the model to be processed, the update amount of the velocity model, and the update amount of the resistivity model, obtain the current velocity model and the current resistivity model; Step S800: Take the current velocity model and the current resistivity model as the models to be processed, and return to the steps of obtaining the target forward modeling response data of the cross-aperture radar travel time CT and the high-density resistivity method based on the models to be processed, until the number of iterations reaches the maximum value, or the fitting difference between the target forward modeling response data and the observation data reaches the preset threshold, and output the joint inversion result of alternating iterations.

[0030] In step S100 of some embodiments, electromagnetic wave data and apparent resistivity data in the target observation area are obtained by pre-deploying high-density resistivity survey lines and receiving and transmitting antennas of cross-hole radar in underground karst, providing a multi-source physical field observation basis for joint inversion and ensuring the completeness and authenticity of the inversion input data.

[0031] In some embodiments, step S100 may include, but is not limited to, steps S110 to S130: Step S110: Obtain apparent resistivity data based on the high-density resistivity survey lines pre-arranged in the target observation area of ​​underground karst. Step S120: Two boreholes are set up in the target observation area, and a cross-hole radar receiving antenna and a cross-hole radar transmitting antenna are set up in the two boreholes respectively to form an overlapping observation area. Step S130: Based on the cross-hole radar receiving antenna and cross-hole radar transmitting antenna in the overlapping observation area of ​​underground karst, electromagnetic wave data is collected by moving along the depth direction of the borehole at a fixed measuring point spacing. Among them, the distance between measuring points is less than half the wavelength of the main frequency of the cross-hole radar.

[0032] In step S110 of some embodiments, apparent resistivity data is obtained based on high-density resistivity survey lines pre-arranged in the underground karst target observation area.

[0033] In step S120 of some embodiments, two boreholes are deployed in the observation area of ​​the underground karst target to be detected, and the measurement depth is determined according to the site conditions. The receiving antenna and transmitting antenna of the cross-hole radar are respectively installed in different boreholes. Figure 2 As shown, the observation area for underground karst targets is the entire high-density resistivity model area, while the area measured by the cross-hole radar CT method is only a part of the entire high-density resistivity model area. Therefore, the area where these two methods (high-density resistivity method and cross-hole radar CT method) overlap in measurement locations is called the overlapping observation area, i.e., the cross-hole radar velocity model. For example, refer to... Figure 2Holes #1 and #2 were laid out in the entire high-density resistivity model region, forming an overlapping observation area for cross-hole radar CT measurement and a non-overlapping observation area in the high-density resistivity model region.

[0034] In step S130 of some embodiments, the receiving antenna and transmitting antenna of the cross-aperture radar are located in different boreholes and move along the borehole depth direction with a fixed measurement point spacing to collect electromagnetic wave data. The measurement point spacing is less than half the wavelength of the cross-aperture radar main frequency.

[0035] In step S200 of some embodiments, an objective function for joint inversion by alternating iteration of cross-hole radar travel time CT and high-density resistivity method with local structural constraints is constructed. By introducing an alternating iteration strategy, the ambiguity of the inversion results of a single method is effectively reduced without relying on complex rock physical property relationships.

[0036] In some embodiments, step S200 may include, but is not limited to, steps S210 to S250: Step S210: Based on the model gradient within the overlapping observation area of ​​the cross-aperture radar travel-time CT and the high-density resistivity method, construct local structural constraints. Step S220: Obtain observation data from the transaperture radar travel-time CT and high-density resistivity method in the electromagnetic wave data and apparent resistivity data; Step S230: Obtain the inverse matrix of the standard deviation of the observation data to obtain the data weighting matrix of the cross-aperture radar travel time CT and the high-density resistivity method; Step S240: Obtain the discretized inversion grid structure of the target observation area. Based on the discretized inversion grid structure, construct a second-order difference operator to obtain the model weighting matrix of the cross-aperture radar travel time CT and the high-density resistivity method. Step S250: Based on the local structural constraints, observation data, data weighting matrix, and model weighting matrix, construct the objective function for the alternating iterative joint inversion of the cross-aperture radar travel-time CT and the high-density resistivity method under local structural constraints.

[0037] In step S210 of some embodiments, the model gradient within the overlapping observation area of ​​the cross-aperture radar travel-time CT is obtained. And obtain the model gradient within the overlapping observation region of the high-density resistivity method. Based on the model gradient within the overlapping observation area of ​​the transapsette radar travel-time CT and the high-density resistivity method, the local structural constraints can be constructed as follows: ; In the formula, The normalized cross gradient function represents the travel time CT of the transap radar and the high-density resistivity method; Represents the current velocity model; This represents the current resistivity model; This represents the velocity model within the overlapping observation area.

[0038] In step S220 of some embodiments, observation data of cross-aperture radar travel-time CT in electromagnetic wave data is acquired. And observational data from the high-density resistivity method within the apparent resistivity data. .

[0039] In step S230 of some embodiments, observation data of the cross-aperture radar travel-time CT is acquired. The reciprocal of the standard deviation matrix is ​​used to obtain the data weighting matrix of the cross-aperture radar travel-time CT. And obtain observation data using the high-density resistivity method. The reciprocal matrix of the standard deviation is used to obtain the data weighting matrix for the high-density resistivity method. .

[0040] In step S240 of some embodiments, the discretized inversion grid structure of the target observation area is obtained. Based on the discretized inversion grid structure, a second-order difference operator (Laplace operator) is constructed to obtain the model weighting matrix of the cross-aperture radar travel-time CT. and the model weighting matrix of the high-density resistivity method .

[0041] In step S250 of some embodiments, according to local structural constraints Observational data and Data weighting matrix and and model weighting matrix and The objective function, constructed using alternating iterations of local structural constraints for transap-hole radar travel-time CT and high-density resistivity method, includes the following formulas:

[0042] ;

[0043] ; In the formula, The objective function representing the travel time CT of the transap-hole radar; The objective function representing the high-density resistivity method; The data fitting term representing the travel-time CT of the transap-hole radar; The data fitting term represents the high-density resistivity method; The structural constraint term representing the travel time CT of the transap-hole radar; The structural constraint term representing the high-density resistivity method; This represents the cross gradient term between the cross-aperture radar travel-time CT and the high-density resistivity method. The adaptive regularization factor representing the travel-time CT of the cross-aperture radar; The adaptive regularization factor representing the high-density resistivity method; The weighting factor of the cross gradient term represents the cross-gradient term of the cross-aperture radar travel-time CT; The weighting factor for the cross gradient term in the high-density resistivity method; This represents the first forward modeling response data of the cross-aperture radar travel-time CT obtained by performing forward modeling on the current velocity model. This represents the second forward modeling response data obtained by performing forward modeling on the current resistivity model using the high-density resistivity method. Represents the initial velocity model; Represents the initial resistivity model; This represents the gradient of the current velocity model; This represents the gradient of the current resistivity model; It represents the square of the L2 norm of the vector.

[0044] In step S300 of some embodiments, the uniform half-space model (high-density resistivity model) is divided into rectangular mesh cells along the x and z coordinate axes in the Cartesian coordinate system, respectively, and the initial velocity model... (like Figure 3 (b) and the initial resistivity model (like Figure 3 (a) is set to the average velocity and average resistivity within the measurement area, respectively.

[0045] In some embodiments, step S300 may include, but is not limited to, steps S310 to S320: Step S310: Obtain the average velocity within the overlapping observation area and set the average velocity as the initial velocity model; Step S320: Obtain the average resistivity within the target observation area and set the average resistivity as the initial resistivity model.

[0046] In step S310 of some embodiments, the target observation area is divided into rectangular grid cells along the x and z coordinate axes in the Cartesian coordinate system, and the initial velocity model is set as the average velocity within the overlapping observation area of ​​the target observation area.

[0047] In step S320 of some embodiments, the target observation area is divided into rectangular grid cells along the x and z coordinate axes in the Cartesian coordinate system, and the initial resistivity model is set as the average resistivity within the target observation area.

[0048] In step S400 of some embodiments, the forward modeling calculation of the cross-aperture radar uses the finite difference method to solve the equation of the process function, ultimately obtaining the first forward response data; the forward modeling calculation of the high-density resistivity method uses the finite element method to calculate the forward response, obtaining the second forward response data; the Jacobian matrix is ​​solved using the interchange theorem. Optionally, the first forward response data... Second forward response data Expanding the series by a first-order Taylor series, we have the following expression: ; ; In the formula, This represents the first forward modeling response data of the cross-aperture radar travel-time CT obtained by performing forward modeling on the initial velocity model. This represents the second forward modeling response data obtained by performing forward modeling on the initial resistivity model using the high-density resistivity method. The first Jacobian matrix representing the travel time CT of the transap radar; This represents the second Jacobian matrix of the high-density resistivity method.

[0049] In some embodiments, step S400 may include, but is not limited to, steps S410 to S420: Step S410: Solve the equation of the process function using the finite difference method, and perform forward modeling response calculation on the travel time CT of the cross-aperture radar to obtain the first forward modeling response data of the travel time CT of the cross-aperture radar. Step S420: The finite element method is used to perform forward response calculation on the high-density resistivity method to obtain the second forward response data of the high-density resistivity method. The target forward response data includes the first forward response data and the second forward response data.

[0050] In step S410 of some embodiments, the first forward response data of the cross-aperture radar travel time CT is obtained by solving the equation of the process function using the finite difference method to calculate the forward response of the cross-aperture radar travel time CT.

[0051] In step S420 of some embodiments, the second forward response data of the high-density resistivity method is obtained by performing forward response calculation of the high-density resistivity method using the finite element method.

[0052] In some embodiments, step S400 may also include, but is not limited to, steps S430 to S440: Step S430: Using the interchange theorem, obtain the first Jacobian matrix of the cross-aperture radar travel time CT; Step S440: Using the interchange theorem, obtain the second Jacobian matrix of the high-density resistivity method; The target Jacobian matrix includes the first Jacobian matrix and the second Jacobian matrix.

[0053] In steps S430 to S440 of some embodiments, the first Jacobian matrix of the cross-aperture radar travel time CT and the second Jacobian matrix of the high-density resistivity method are both obtained by solving using the interchange theorem.

[0054] In step S500 of some embodiments, the normalized cross gradient partial derivatives within the overlapping observation region can be calculated by minimizing the objective function, as expressed below: ; ; ; ; In the formula, The global cross gradient partial derivative represents the cross-aperture radar travel time CT. The global cross gradient partial derivative represents the high-density resistivity method; The partial derivative of the local cross gradient represents the cross-aperture CT of the cross-aperture radar travel time. The partial derivatives of the local cross gradients in the high-density resistivity method; Represents the velocity model within the overlapping observation region; Resistivity model representing the overlapping observation region; Normalized cross gradient function representing the travel time CT of transap radar and the high-density resistivity method .

[0055] In step S600 of some embodiments, a large matrix of cross-aperture radar travel time CT and high-density resistivity method is established, and the least squares method is used to solve the velocity model update and resistivity model update of alternating iterative joint inversion respectively.

[0056] In some embodiments, step S600 may include, but is not limited to, steps S610 to S650: Step S610: Obtain observation data from the transaperture radar travel-time CT and high-density resistivity method in the electromagnetic wave data and apparent resistivity data; Step S620: Obtain the inverse matrix of the standard deviation of the observation data to obtain the data weighting matrix of the cross-aperture radar travel time CT and the high-density resistivity method; Step S630: Obtain the residual between the observed data and the target forward response data; Step S640: Obtain the discretized inversion grid structure of the target observation area. Based on the discretized inversion grid structure, construct a second-order difference operator to obtain the model weighting matrix of the cross-aperture radar travel time CT and the high-density resistivity method. Step S650: Construct a large matrix based on the data weighting matrix, residuals, model weighting matrix, target Jacobian matrix, and cross-gradient partial derivatives.

[0057] In step S610 of some embodiments, observation data of the transap-hole radar travel-time CT in the electromagnetic wave data is acquired. And observational data from the high-density resistivity method within the apparent resistivity data. .

[0058] In step S620 of some embodiments, observation data of the cross-aperture radar travel-time CT is acquired. The reciprocal of the standard deviation matrix is ​​used to obtain the data weighting matrix of the cross-aperture radar travel-time CT. And obtain observation data using the high-density resistivity method. The reciprocal matrix of the standard deviation is used to obtain the data weighting matrix for the high-density resistivity method. .

[0059] In step S630 of some embodiments, the residual between the observed data and the target forward response data is calculated, using a formula including: ; ; In the formula, The residual between the observation data of the cross-aperture radar travel-time CT and the first forward response data; The residual represents the difference between the observation data from the high-density resistivity method and the second forward response data.

[0060] In step S640 of some embodiments, a second-order difference operator (Laplace operator) is constructed based on the discretized inversion grid structure of the target observation area, which yields the model weighting matrix of the cross-aperture radar travel-time CT. and the model weighting matrix of the high-density resistivity method .

[0061] In step S650 of some embodiments, a large matrix can be constructed based on the data weighting matrix, residuals, model weighting matrix, target Jacobian matrix, and cross-gradient partial derivatives, as expressed below: ; ; In the formula, The rate of change of the model represents the rate of change of alternating iterative joint inversion; This represents the update amount of the resistivity model in the alternating iterative joint inversion.

[0062] In step S700 of some embodiments, the current velocity model and the current resistivity model are obtained based on the model to be processed, the velocity model update amount, and the resistivity model update amount. For example, the current velocity model can be obtained by adding the velocity model update amount to the initial velocity model or the velocity model obtained in the previous iteration; the current resistivity model can be obtained by adding the resistivity model update amount to the resistivity model or the resistivity model obtained in the previous iteration. Taking the initial velocity model and the initial resistivity model as the models to be processed as an example, the formulas for calculating the current velocity model and the current resistivity model obtained from alternating iterative joint inversion include: ; .

[0063] In step S800 of some embodiments, the current velocity model and current resistivity model obtained in the current iteration are used as models to be processed. The process returns to the step of obtaining the target forward response data using the cross-aperture radar travel-time CT and the high-density resistivity method based on the models to be processed, and calculating the fitting difference between the target forward response data and the observed data. If the number of iterations reaches a maximum value or the fitting difference reaches a threshold, the iterative inversion stops, and the final alternating iterative joint inversion result is output. For example, the formula for calculating the fitting difference between the target forward response data and the observed data includes: ; ; In the formula, This represents the fitting difference between the first forward response data and the observation data from the cross-aperture radar travel-time CT. This represents the fitting difference between the second forward response data and the observation data from the high-density resistivity method; Represents the total number of observation data from the transap-hole radar travel-time CT. This is an index of observation data from cross-aperture radar travel-time CT. This represents the total number of observations using the high-density resistivity method. This is an index of observational data for the high-density resistivity method.

[0064] like Figure 4 As shown, the specific implementation steps of the joint inversion of cross-aperture radar travel time CT and high-density resistivity alternating iteration include: a. Acquire electromagnetic wave data and apparent resistivity data within the observation area; b. Construct the objective function by alternating iterative joint inversion of cross-aperture radar travel time CT and high-density resistivity with local structural constraints; c. Establish the initial velocity model and the initial resistivity model; d. Perform forward modeling response calculations using the cross-aperture radar travel time CT and high-density resistivity method, and solve for the Jacobian matrix on the initial or iterative model; e. Calculate the normalized cross gradient partial derivatives within the overlapping observation region; f. Establish a large matrix for the cross-aperture radar travel time CT and the high-density resistivity method, and use the least squares method to solve for the velocity and resistivity updates obtained by alternating iterative joint inversion; g. Calculate the velocity and resistivity iterative model for alternating iterative joint inversion; h) Perform forward modeling calculations on the velocity and resistivity iterative models respectively, and calculate the fitting difference between the forward response data and the observed data. If the number of iterations reaches the maximum value or the fitting difference reaches the threshold, stop the iterative inversion and output the final alternating iterative joint inversion result; otherwise, return to step d.

[0065] In some embodiments, the accuracy and effectiveness of the joint inversion method of alternating iterative retrieval of cross-aperture radar travel time CT and high-density resistivity under local structural constraints are verified. This embodiment designs as follows: Figure 3 The low resistance shown The low-velocity anomaly model performs individual and joint inversions on the synthetic data, and compares and analyzes the velocity and resistivity inversion results.

[0066] The resistivity model in this embodiment has geometric dimensions of 160m × 36m; the mesh element size is 1m × 1m; the resistivity of the background medium is 900 Ω·m; a circular low-resistivity anomaly with a resistivity of 100 Ω·m and a diameter of 8m is set in the background medium, with its center located at 80m and 18m; in the resistivity forward modeling, the acquisition device adopts a dipole-dipole electrode array with a unit electrode spacing of 2m, containing a total of 81 electrodes. The model extension region set to suppress boundary effects is not included in the model extension region. Figure 3 As shown in section (a). The two boreholes of the velocity model in the embodiment are located at 72m and 88m along the x-axis of the resistivity model, respectively. The velocity model has a geometry of 16m × 34m; the mesh cell size is 0.05m × 0.05m; the relative permittivity of the background medium is 11; and the relative permittivity of the circular low-velocity anomaly is 25. Figure 3Part (b) also does not show the PML layer set by the absorbing boundary conditions. In the velocity forward modeling simulation, the transmitting and receiving antennas are respectively placed in two boreholes, which are marked by circles and triangles in the schematic diagram. The transmitting antenna excites the electromagnetic wave source at a fixed depth, while the receiving antenna moves in a step size of 1m in another borehole to achieve dense signal acquisition with one transmission and multiple receptions. The center frequency of the Ricker wavelet is set to 50MHz, the sampling interval is 0.118ns, and the recording length is 6000 sampling points. Whether it is a single inversion or a joint inversion, the maximum number of iterations is set to 30, and the initial model is a uniform half-space model. The inversion grid sizes of the velocity and resistivity models are 0.5m×0.5m and 1m×1m, respectively. The structural similarity index (SSIM), fitting error (RMS), and average resistivity (ρ) will be used to calculate the structural similarity index (SSIM), fitting error (RMS), and average resistivity (ρ). mean ) and average velocity (V) mean The inversion imaging effect is evaluated from four aspects.

[0067] When all cross gradient weight factors are set to 0, the corresponding Figure 5 The individual inversion results; among them, Figure 5 Parts (a), (b), (c), and (d) in the diagram represent the velocity inversion result, the resistivity inversion result, the cross-gradient plot, and the fitting error convergence curve, respectively. When the cross-gradient weight factor is set to 20-10, the corresponding... Figure 6 Alternating iterative joint inversion; where, Figure 6 Parts (a), (b), (c), and (d) in the figure represent the joint velocity inversion results, the joint resistivity inversion results, the cross gradient plot, and the fitting error convergence curve, respectively. Regardless of whether the inversion method is used individually or jointly, both cross-aperture radar travel-time CT and the high-density resistivity method can reliably identify the horizontal center position of anomalies, but they differ significantly in structural morphology reconstruction and physical property parameter recovery.

[0068] like Figure 5 As shown in part (a), the travel time results of the transaperture radar retrieved separately produce obvious trailing artifacts in the velocity model due to insufficient ray coverage, resulting in a low structural similarity index (SSIM). CT =0.781). However, the average velocity value (V) of the reconstructed low-velocity anomaly mean =0.077m / ns) is still relatively close to the actual model. For example Figure 5 As shown in section (b), the average resistivity ρ of the high-density resistivity method retrieved alone. mean Structural Similarity Index (SSIM) ERT The values ​​are 352.99 Ω·m and 0.572, respectively, but the anomaly exhibits a false anomaly phenomenon extending downwards, making it difficult to accurately characterize its longitudinal boundary. Furthermore, the cross-gradient values ​​between the two models have a wide distribution range (-6 × 10⁻⁶).- ³ to +6×10 - ³ indicates that the structural consistency between individually inverted models is poor.

[0069] like Figure 6 As shown in part (a), the background resistivity structure information provided by the resistivity model can effectively suppress the trailing artifacts caused by insufficient ray coverage during the travel time of the transaperture radar, making the morphology of low-velocity anomalies closer to the true velocity model. Its structural similarity index SSIM CT Increased to 0.896. For example... Figure 6 As shown in part (b), the joint inversion effectively utilizes high-resolution transap-scan radar travel time structure information to accurately characterize the longitudinal top and bottom boundaries of the low resistivity anomaly, making the anomaly's shape and location closer to the true resistivity model. Although the average resistivity (ρ) of the resistivity model... mean =355.06Ω·m) and structural similarity index (SSIM) ERT =0.569) slightly decreased. The cross gradient values ​​between the two models range from -2×10 - ³ to +2×10 - The results are between 3 and 4, significantly lower than the results of the individual inversion, which verifies that the joint inversion achieves better results in both spatial morphology and numerical recovery of physical properties.

[0070] This invention also provides a joint inversion device for alternating iterative inversion of cross-aperture radar travel time CT and high-density resistivity, which can realize the above-mentioned joint inversion method for alternating iterative inversion of cross-aperture radar travel time CT and high-density resistivity. The device includes: The data acquisition module is used to acquire electromagnetic wave data and apparent resistivity data within the target observation area; The objective function construction module is used to construct an objective function that is jointly inverted by alternating iterations of cross-aperture radar travel-time CT and high-density resistivity method based on electromagnetic wave data and apparent resistivity data. The initial model building module is used to establish the initial velocity model and the initial resistivity model; the initial velocity model and the initial resistivity model are used as models to be processed. The forward modeling module is used to obtain the target forward modeling response data of the cross-aperture radar travel time CT and the high-density resistivity method based on the model to be processed, and to obtain the target Jacobian matrix of the cross-aperture radar travel time CT and the high-density resistivity method based on the target forward modeling response data. The cross-gradient partial derivative acquisition module is used to obtain the cross-gradient partial derivatives of the overlapping observation regions within the target observation region by minimizing the objective function. The model update acquisition module is used to establish a large matrix for cross-aperture radar travel time CT and high-density resistivity method based on target forward response data, target Jacobian matrix and cross gradient partial derivatives, and to acquire the velocity model update and resistivity model update based on the large matrix. The model update module is used to obtain the current velocity model and the current resistivity model based on the model to be processed, the update amount of the velocity model, and the update amount of the resistivity model. The iterative inversion module is used to take the current velocity model and the current resistivity model as the models to be processed, and return the steps of obtaining the target forward modeling response data of the cross-aperture radar travel time CT and the high-density resistivity method based on the models to be processed, until the number of iterations reaches the maximum value, or the fitting difference between the target forward modeling response data and the observation data reaches the preset threshold, and outputs the joint inversion results of alternating iterations.

[0071] It is understood that the content of the above method embodiments is applicable to the present device embodiments. The specific functions implemented by the present device embodiments are the same as those of the above method embodiments, and the beneficial effects achieved are also the same as those achieved by the above method embodiments.

[0072] This invention also provides an electronic device, which includes a processor and a memory. The memory stores a computer program, and the processor executes the computer program to implement the above-described method. This electronic device can be any smart terminal, including a tablet computer, an in-vehicle computer, or similar device.

[0073] It is understood that the content of the above method embodiments is applicable to this device embodiment. The specific functions implemented by this device embodiment are the same as those of the above method embodiments, and the beneficial effects achieved are also the same as those achieved by the above method embodiments.

[0074] refer to Figure 7 , Figure 7 The hardware structure of an electronic device according to another embodiment is illustrated. The electronic device includes: The processor 901 can be implemented using a general-purpose central processing unit (CPU), microprocessor, application specific integrated circuit (ASIC), or one or more integrated circuits, and is used to execute relevant programs to implement the technical solutions provided in the embodiments of the present invention. The memory 902 can be implemented as a read-only memory (ROM), static storage device, dynamic storage device, or random access memory (RAM). The memory 902 can store the operating system and other application programs. When the technical solutions provided in the embodiments of this specification are implemented through software or firmware, the relevant program code is stored in the memory 902 and is called and executed by the processor 901. The input / output interface 903 is used to implement information input and output; The communication interface 904 is used to enable communication and interaction between this device and other devices. Communication can be achieved through wired means (such as USB, Ethernet cable, etc.) or wireless means (such as mobile network, WIFI, Bluetooth, etc.). Bus 905 transmits information between various components of the device (e.g., processor 901, memory 902, input / output interface 903, and communication interface 904); The processor 901, memory 902, input / output interface 903, and communication interface 904 are connected to each other within the device via bus 905.

[0075] This invention also provides a computer-readable storage medium storing a computer program that, when executed by a processor, implements the above-described method.

[0076] It is understood that the content of the above method embodiments is applicable to this storage medium embodiment. The specific functions implemented in this storage medium embodiment are the same as those in the above method embodiments, and the beneficial effects achieved are also the same as those achieved in the above method embodiments.

[0077] This invention also provides a computer program product or computer program that includes computer instructions stored in a computer-readable storage medium. A processor of a computer device can read the computer instructions from the computer-readable storage medium and execute the computer instructions to cause the computer device to perform the aforementioned method.

[0078] In summary, the cross-aperture radar travel-time CT and high-density resistivity alternating iterative joint inversion method and related equipment of the present invention have the following advantages: 1. This invention introduces an alternating iterative joint inversion strategy into a joint inversion method of cross-hole radar travel-time CT and high-density resistivity under local structural constraints. Compared to individual inversions, this method does not rely on complex rock property relationships and can effectively reduce the ambiguity of the inversion results.

[0079] 2. In this embodiment of the invention, an adaptive regularization factor is set for the structural constraint term to balance the data fitting term and the structural constraint term; a cross gradient weight factor is set for the cross gradient term to balance the constraint effect of the velocity and resistivity iterative model, thereby improving the reliability of the inversion results.

[0080] 3. The embodiments of the present invention can effectively avoid inversion errors caused by direct coupling of physical property parameters with large differences in magnitude by normalizing the cross gradient function and partial derivatives.

[0081] In some alternative embodiments, the functions / operations mentioned in the block diagrams may not occur in the order shown in the operation diagrams. For example, depending on the functions / operations involved, two consecutively shown blocks may actually be executed substantially simultaneously, or the blocks may sometimes be executed in reverse order. Furthermore, the embodiments presented and described in the flowcharts of this invention are provided by way of example to provide a more comprehensive understanding of the technology. The disclosed methods are not limited to the operations and logic flows presented herein. Alternative embodiments are contemplated in which the order of various operations is altered and sub-operations described as part of a larger operation are executed independently.

[0082] Furthermore, although the invention has been described in the context of functional modules, it should be understood that, unless otherwise stated, one or more of the described functions and / or features may be integrated into a single physical device and / or software module, or one or more functions and / or features may be implemented in a separate physical device or software module. It is also understood that a detailed discussion of the actual implementation of each module is unnecessary for understanding the invention. Rather, given the properties, functions, and internal relationships of the various functional modules in the apparatus disclosed herein, the actual implementation of the module will be understood within the scope of conventional skill of an engineer. Therefore, those skilled in the art can implement the invention as set forth in the claims using ordinary techniques without excessive experimentation. It is also understood that the specific concepts disclosed are merely illustrative and not intended to limit the scope of the invention, which is determined by the full scope of the appended claims and their equivalents.

[0083] If the aforementioned functions are implemented as software functional units and sold or used as independent products, they can be stored in a computer-readable storage medium. Based on this understanding, the technical solution of this invention, or the part that contributes to the prior art, or a portion of the technical solution, can be embodied in the form of a software product. This computer software product is stored in a storage medium and includes several instructions to cause a computer device (which may be a personal computer, server, or network device, etc.) to execute all or part of the steps of the methods described in the various embodiments of this invention. The aforementioned storage medium includes various media capable of storing program code, such as USB flash drives, portable hard drives, read-only memory, random access memory, magnetic disks, or optical disks.

[0084] The logic and / or steps represented in the flowchart or otherwise described herein, for example, can be considered as a sequenced list of executable instructions for implementing logical functions, and can be embodied in any computer-readable medium for use by, or in conjunction with, an instruction execution system, apparatus, or device (such as a computer-based system, a processor-including system, or other system that can fetch and execute instructions from, an instruction execution system, apparatus, or device). For the purposes of this specification, "computer-readable medium" can be any means that can contain, store, communicate, propagate, or transmit programs for use by, or in conjunction with, an instruction execution system, apparatus, or device.

[0085] It should be understood that various parts of the present invention can be implemented in hardware, software, firmware, or a combination thereof. In the above embodiments, multiple steps or methods can be implemented in software or firmware stored in memory and executed by a suitable instruction execution system. For example, if implemented in hardware, as in another embodiment, it can be implemented using any one or a combination of the following techniques known in the art: discrete logic circuits having logic gates for implementing logical functions on data signals, application-specific integrated circuits (ASICs) having suitable combinational logic gates, programmable gate arrays (PGAs), field-programmable gate arrays (FPGAs), etc.

[0086] In the description of this specification, references to terms such as "one embodiment," "some embodiments," "example," "specific example," or "some examples," etc., indicate that a specific feature, structure, material, or characteristic described in connection with that embodiment or example is included in at least one embodiment or example of the invention. In this specification, the illustrative expressions of the above terms do not necessarily refer to the same embodiment or example. Furthermore, the specific features, structures, materials, or characteristics described may be combined in any suitable manner in one or more embodiments or examples.

[0087] Although embodiments of the invention have been shown and described, those skilled in the art will understand that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the claims and their equivalents.

[0088] The above is a detailed description of the preferred embodiments of the present invention. However, the present invention is not limited to the embodiments described. Those skilled in the art can make various equivalent modifications or substitutions without departing from the spirit of the present invention. All such equivalent modifications or substitutions are included within the scope defined by the claims of the present invention.

Claims

1. A joint inversion method for alternating iterative inversion of cross-aperture radar travel time CT and high-density resistivity, characterized in that, Includes the following steps: Acquire electromagnetic wave data and apparent resistivity data within the target observation area; Based on the electromagnetic wave data and the apparent resistivity data, an objective function is constructed by alternating iterative joint inversion of cross-aperture radar travel-time CT and high-density resistivity method. Establish an initial velocity model and an initial resistivity model; use the initial velocity model and the initial resistivity model as models to be processed; Based on the model to be processed, obtain the target forward modeling response data of the cross-aperture radar travel time CT and the high-density resistivity method, and based on the target forward modeling response data, obtain the target Jacobian matrix of the cross-aperture radar travel time CT and the high-density resistivity method. By minimizing the objective function, the cross gradient partial derivatives of the overlapping observation regions within the target observation region are obtained; Based on the target forward response data, the target Jacobian matrix, and the cross gradient partial derivatives, a large matrix is ​​established for the cross-aperture radar travel time CT and the high-density resistivity method, and the velocity model update and resistivity model update are obtained based on the large matrix. Based on the model to be processed, the velocity model update amount, and the resistivity model update amount, obtain the current velocity model and the current resistivity model; The current velocity model and the current resistivity model are used as the model to be processed. The process of obtaining the target forward modeling response data of the cross-aperture radar travel time CT and the high-density resistivity method based on the model to be processed is repeated until the number of iterations reaches the maximum value or the fitting difference between the target forward modeling response data and the observation data reaches a preset threshold. The joint inversion result of alternating iterations is then output.

2. The method according to claim 1, characterized in that, The acquisition of electromagnetic wave data and apparent resistivity data within the target observation area includes the following steps: The apparent resistivity data is obtained based on the high-density resistivity survey lines pre-arranged in the target observation area of ​​the underground karst. Two boreholes are drilled in the target observation area, and a cross-hole radar receiving antenna and a cross-hole radar transmitting antenna are respectively installed in the two boreholes to form the overlapping observation area; Based on the cross-hole radar receiving antenna and the cross-hole radar transmitting antenna in the overlapping observation area of ​​the underground karst, the electromagnetic wave data is collected by moving along the depth direction of the borehole at a fixed measuring point spacing. The distance between the measuring points is less than half the wavelength of the main frequency of the cross-hole radar.

3. The method according to claim 1, characterized in that, The objective function constructed based on the electromagnetic wave data and the apparent resistivity data, through alternating iterative joint inversion of cross-aperture radar travel-time CT and high-density resistivity method, includes the following steps: Based on the model gradient within the overlapping observation area using cross-aperture radar travel-time CT and high-density resistivity method, local structural constraints are constructed. Acquire the observation data from the electromagnetic wave data and the apparent resistivity data, including the cross-aperture radar travel-time CT and the high-density resistivity method. Obtain the reciprocal matrix of the standard deviation of the observation data to get the data weighting matrix of cross-aperture radar travel time CT and high-density resistivity method; The discretized inversion grid structure of the target observation area is obtained. Based on the discretized inversion grid structure, a second-order difference operator is constructed to obtain the model weighting matrix of the cross-aperture radar travel time CT and the high-density resistivity method. Based on the local structural constraints, the observation data, the data weighting matrix, and the model weighting matrix, the objective function is constructed by alternating iterations of the cross-aperture radar travel-time CT and the high-density resistivity method under local structural constraints.

4. The method according to claim 1, characterized in that, The establishment of the initial velocity model and the initial resistivity model includes the following steps: Obtain the average velocity within the overlapping observation area, and set the average velocity as the initial velocity model; Obtain the average resistivity within the target observation area, and set the average resistivity as the initial resistivity model.

5. The method according to claim 1, characterized in that, The process of obtaining target forward modeling response data based on the model to be processed, using cross-aperture radar travel-time CT and high-density resistivity method, includes the following steps: The finite difference method is used to solve the equation of the process function, and the forward modeling response of the cross-aperture radar travel-time CT is calculated to obtain the first forward modeling response data of the cross-aperture radar travel-time CT. The finite element method is used to perform forward response calculations on the high-density resistivity method to obtain the second forward response data of the high-density resistivity method; The target forward modeling response data includes the first forward modeling response data and the second forward modeling response data.

6. The method according to claim 1, characterized in that, The step of obtaining the target Jacobian matrix using the cross-aperture radar travel-time CT and the high-density resistivity method based on the target forward response data includes the following steps: The first Jacobian matrix of the cross-aperture radar travel time CT is obtained by using the interchange theorem; The second Jacobian matrix of the high-density resistivity method is obtained by using the interchange theorem. The target Jacobian matrix includes the first Jacobian matrix and the second Jacobian matrix.

7. The method according to claim 1, characterized in that, The process of establishing a large matrix for cross-aperture radar travel-time CT and high-density resistivity method based on the target forward response data, the target Jacobian matrix, and the cross gradient partial derivatives includes the following steps: Acquire the observation data from the electromagnetic wave data and the apparent resistivity data, including the cross-aperture radar travel-time CT and the high-density resistivity method. Obtain the reciprocal matrix of the standard deviation of the observation data to get the data weighting matrix of cross-aperture radar travel time CT and high-density resistivity method; Obtain the residual between the observed data and the target forward response data; The discretized inversion grid structure of the target observation area is obtained. Based on the discretized inversion grid structure, a second-order difference operator is constructed to obtain the model weighting matrix of the cross-aperture radar travel time CT and the high-density resistivity method. The large matrix is ​​constructed based on the data weighting matrix, the residual, the model weighting matrix, the target Jacobian matrix, and the cross gradient partial derivatives.

8. A joint inversion device for alternating iterative inversion of cross-aperture radar travel time CT and high-density resistivity, characterized in that, include: The data acquisition module is used to acquire electromagnetic wave data and apparent resistivity data within the target observation area; The objective function construction module is used to construct an objective function jointly inverted by alternating iterations of cross-aperture radar travel-time CT and high-density resistivity method based on the electromagnetic wave data and the apparent resistivity data. An initial model building module is used to establish an initial velocity model and an initial resistivity model; the initial velocity model and the initial resistivity model are used as models to be processed. The forward modeling module is used to obtain the target forward modeling response data of the cross-aperture radar travel time CT and the high-density resistivity method based on the model to be processed, and to obtain the target Jacobian matrix of the cross-aperture radar travel time CT and the high-density resistivity method based on the target forward modeling response data. The cross-gradient partial derivative acquisition module is used to obtain the cross-gradient partial derivatives of the overlapping observation regions within the target observation region by minimizing the objective function; The model update quantity acquisition module is used to establish a large matrix of cross-aperture radar travel time CT and high-density resistivity method based on the target forward response data, the target Jacobian matrix and the cross gradient partial derivatives, and to acquire the velocity model update quantity and resistivity model update quantity based on the large matrix. The model update module is used to obtain the current velocity model and the current resistivity model based on the model to be processed, the velocity model update amount, and the resistivity model update amount. The iterative inversion module is used to take the current velocity model and the current resistivity model as the model to be processed, and return the step of obtaining the target forward modeling response data of the cross-aperture radar travel time CT and the high-density resistivity method according to the model to be processed, until the number of iterations reaches the maximum value, or the fitting difference between the target forward modeling response data and the observation data reaches a preset threshold, and output the joint inversion result of alternating iterations.

9. An electronic device, characterized in that, Including the processor and memory; The memory is used to store programs; The processor executes the program to implement the method as described in any one of claims 1 to 7.

10. A computer program product, comprising a computer program, characterized in that, When the computer program is executed by a processor, it implements the method as described in any one of claims 1 to 7.

Images