A frequency domain marine controlled source electromagnetic and magnetotelluric data joint imaging method

Abstract

This invention discloses a joint imaging method for frequency-domain controlled-source electromagnetic (CMY) and magnetotelluric (MT) data from the ocean. Its key feature is that, addressing the shortcomings of conventional electromagnetic data imaging methods such as insufficient imaging accuracy and low resolution, this invention first employs an unstructured dual-grid system to partition the model, using one set as a forward grid and the other as an inversion grid. Then, the Fast Occam Imaging Algorithm is used to achieve 3D inversion imaging of CMY and MT data separately. Finally, the CMY and MT data are added to the same inversion dataset, and a normalized joint fitting function is introduced to ensure that the inversion results can evenly fit different subsets of electromagnetic data, resulting in a new 3D joint imaging method for CMY and MT data. This invention overcomes the shortcomings of inversion using a single electromagnetic dataset, improves the accuracy and resolution of electromagnetic data inversion imaging, and provides users with a powerful tool.

Description

Technical Field

[0001] This invention relates to the field of combined imaging methods using multiple geophysical data, specifically a combined imaging method using frequency-domain ocean controlled-source electromagnetic and magnetotelluric data. Technical Background

[0002] With the continuous deepening of my country's industrialization and urbanization, and its ongoing development towards becoming a modern and powerful nation, the country's overall dependence on foreign mineral resources continues to rise. This is particularly true for oil and gas resources. In 2019, China imported 540 million tons of oil, raising its dependence on foreign sources to 73.5%. In 2020, China imported 101.66 million tons of natural gas, raising its dependence on foreign sources to approximately 84.3%. This high dependence on imported mineral resources has become one of the main bottlenecks restricting my country's development, making the acceleration of the exploration and development of oil, gas, and other mineral resources a crucial national strategic goal. Electromagnetic methods, as a traditional geophysical exploration technique, play a vital role in mineral resource exploration, especially the frequency-domain controlled-source electromagnetic method (FCR). The frequency-domain controlled-source electromagnetic method uses an artificial field source to control the direction, strength, and frequency range of the observed electromagnetic field. It has advantages such as large detection depth, high resolution, and strong anti-interference capabilities in the field, and is widely used in various fields such as mineral resource exploration and environmental engineering exploration. In particular, it has played an important role in the exploration of marine oil and gas resources in recent years. Due to the high sensitivity of frequency domain controllable source electromagnetic method to high resistivity oil and gas reservoirs, it provides a lot of important information for the exploration of marine oil and gas resources.

[0003] Magnetotelluric (MT) is a geophysical method that uses a natural alternating electromagnetic field incident vertically on the Earth's surface at high altitudes as the field source and the electrical differences in rocks as the material basis to study the electrical structure of sedimentary caprock and even the crust and upper mantle. It has advantages such as low cost, no shielding by high-resistivity layers, and high resolution for low-resistivity layers. It has been successfully applied in domestic and international oil and gas exploration and prospecting, the study of the electrical structure of the crust and upper mantle, the investigation of geothermal fields, and the prediction and forecasting of natural earthquakes (Ye Yixin et al., 2020).

[0004] Although there are many imaging methods for controlled-source electromagnetic data and magnetotelluric data, most of them use structured grids to divide the model. For simple geoelectric models, structured grids can achieve high computational accuracy, but for complex geological structures, such as arbitrarily undulating terrain and inclined interfaces, structured grids are difficult to accurately simulate. Moreover, the division and densification of structured grids require users to have rich experience, which is not conducive to imaging complex structures. In addition, existing electromagnetic data imaging methods only consider a single type of electromagnetic data, which has the following defects: (1) there are problems such as insufficient accuracy and low resolution of inversion imaging results; (2) the understanding of complex seabed geological conditions is not deep, and most of the analysis is based on simple models, making it difficult to accurately restore the real seabed structure. Summary of the Invention

[0005] The purpose of this invention is to improve the imaging accuracy and resolution of marine controlled-source electromagnetic and magnetotelluric data imaging methods, and to make up for the shortcomings of single dataset inversion, so as to make the imaging results clearer and higher resolution. A frequency domain marine controlled-source electromagnetic and magnetotelluric data joint imaging method is proposed.

[0006] The frequency-domain joint imaging method for ocean controlled-source electromagnetic (OME) and magnetotelluric (MT) data proposed in this invention is based on the following idea: Addressing the issues of insufficient imaging accuracy and low resolution in conventional OME and MT imaging methods, this invention first uses an unstructured dual-grid system to discretize the forward and inverse model regions, with one set as the forward grid and the other as the inverse grid. Then, the fast Occam's algorithm is used to achieve three-dimensional inverse imaging of OME and MT data respectively. Finally, the OME and MT data are added to the same inverse dataset, and a normalized joint fitting function is introduced to ensure that the inverse results can evenly fit different subsets of electromagnetic data, resulting in a new joint imaging method for OME and MT data.

[0007] The unstructured dual-grid system consists of two sets: one is the inversion grid, i.e., the coarse grid, used for discretizing the parameters of the inversion model; the other is the forward grid, i.e., the fine grid, used for discretizing the forward model to ensure the accuracy of the forward modeling, and is derived from the refinement of the inversion grid.

[0008] The refinement strategy for the forward mesh in this invention is as follows:

[0009] (1) Refine the cells near the receiving point, find the cell where the receiving point is located, and determine whether its volume is less than the set requirement for the volume of the receiving point cell. If it does not meet the requirement, halve its volume and refine it in the subsequent mesh refinement. Repeat the above steps until the volume of the cell meets the requirement.

[0010] (2) For electromagnetic data from a controlled marine source, in order to overcome the singularity of the emission source, the cells near the emission point need to be refined. Locate the cell where the emission point is located and determine whether its volume is smaller than the set requirement for the volume of the emission point cell. If it does not meet the requirement, its volume is halved and refined in the subsequent mesh refinement. Repeat the above steps until the volume of the cell meets the requirement.

[0011] (3) Iterate through each cell of the mesh to find the minimum distance from the cell to the receiver and the transmitter. Based on the distance and skin depth, calculate the maximum volume that the cell should meet. If the volume exceeds the maximum volume, its volume is halved and refined in the subsequent mesh refinement. Repeat the above steps until the volume of all cells meets the requirements.

[0012] (4) For inversion, in order to achieve unbiased transfer of inversion parameters to forward parameters, it is only necessary to mark each unit that needs to be inverted in the inversion grid as an independent region.

[0013] Furthermore, the inversion objective function is constructed using the Tikhonov regularization method, and the functional U is expressed as:

[0014] U=μ||Rm|| 2 +||W(dF(m))|| 2

[0015] In the formula: m is an N-dimensional model parameter vector, typically the resistivity value; R is the roughness operator matrix; μ is the Lagrange multiplier, used to balance model roughness and data fitting error. When μ takes a large value, the inversion mainly searches for a smooth model; otherwise, it mainly searches for the minimum fitting error. W is the associated diagonal weighted matrix; d is the observed data vector; and F(m) is the forward response corresponding to model m.

[0016] Given an initial model m k The objective function is minimized using the following iterative method:

[0017]

[0018] The data correction vector is:

[0019]

[0020] The Jacobian matrix J is an M×N matrix, where each component is expressed as:

[0021]

[0022] Furthermore, the Jacobian matrix is ​​calculated using the adjoint reciprocity theorem:

[0023] Introducing accompanying electrical sources and accompanying magnetic source The following Maxwell's equations are obtained:

[0024]

[0025]

[0026] Considering the boundary conditions and substituting them into the above equation, we obtain the following equation:

[0027]

[0028] The above formula is the basic formula for solving the Jacobian matrix of the electromagnetic field. The Jacobian matrix is ​​calculated by selecting a suitable adjoint field source. If it is necessary to obtain the E at the observation point r0... x Jacobi style, making and Then we get:

[0029]

[0030] The above formula is E x The formula for calculating the Jacobian of a component is used to obtain... The electric field value E can be obtained by first solving for the original electromagnetic field, and then by calculating the electric field value E under the condition of the accompanying source. * Finally, the dot product E * E is sufficient.

[0031] Furthermore, the fitting function in the objective function is extended to:

[0032] ||W(dF(m))|| 2 =||W1(d1-F1(m))|| 2 +||W2(d2-F2(m))|| 2

[0033]

[0034] Assuming the desired fitting error is achieved by increasing the data weighting coefficients α1 and α2, where the values ​​of α1 and α2 are related to the amount of data in each subset, the normalized joint fitting function is as follows:

[0035]

[0036] In the formula It is a data fitting balance weight that normalizes the fitting function for each data subset so that the smaller data subset has the same impact on the overall error as the larger subset.

[0037] Based on the above-mentioned inventive concept, the technical solution proposed by this invention to achieve the invention's objective is a joint imaging method for frequency-domain ocean controllable source electromagnetic and magnetotelluric data, the steps of which are:

[0038] Step 1: Mesh Refinement

[0039] An unstructured dual-grid system is used to discretize the forward and inverse model regions: one set is the inverse grid, i.e., the coarse grid, used to discretize the parameters of the inverse model; the other set is the forward grid, i.e., the fine grid, used to discretize the forward model to ensure the accuracy of the forward model, and is derived from the inverse grid.

[0040] The refinement steps of the forward mesh are as follows:

[0041] (1) Refine the cells near the receiving point, find the cell where the receiving point is located, and determine whether its volume is less than the set requirement for the volume of the receiving point cell. If it does not meet the requirement, its volume is halved and then refined in the subsequent mesh refinement. Repeat the above steps until the volume of the cell meets the requirement.

[0042] (2) If it is ocean-controlled source electromagnetic data, in order to overcome the singularity of the emission source point, the cells near the emission point need to be refined; find the cell where the emission point is located, and determine whether its volume is less than the set requirement for the volume of the emission point cell. If it does not meet the requirement, its volume is halved and refined in the subsequent mesh refinement; repeat the above steps until the volume of the cell meets the requirement.

[0043] (3) Iterate through each cell of the grid to find the minimum distance from the cell to the receiving point and the transmitting source. Based on the distance and skin depth, calculate the maximum volume that the cell should meet. If the volume exceeds the maximum volume, the volume is halved and refined in the subsequent mesh refinement. Repeat the above steps until the volume of all cells meets the requirements.

[0044] (4) For inversion, in order to achieve unbiased transfer of inversion parameters to forward parameters, it is only necessary to mark each unit that needs to be inverted in the inversion grid as an independent region.

[0045] Step 2: Fast Occam Inversion Imaging

[0046] The inversion objective function is constructed using the Tikhonov regularization method, and the functional U is expressed as:

[0047] U=μ||Rm|| 2 +||W(dF(m))|| 2

[0048] In the formula: m is an N-dimensional model parameter vector, typically the resistivity value; R is the roughness operator matrix; μ is the Lagrange multiplier, used to balance model roughness and data fitting error. When μ takes a large value, the inversion mainly searches for a smooth model; otherwise, it mainly searches for the minimum fitting error. W is the associated diagonal weighted matrix; d is the observed data vector; and F(m) is the forward response corresponding to model m.

[0049] Given an initial model m k The objective function is minimized using the following iterative method:

[0050]

[0051] The data correction vector is:

[0052]

[0053] The Jacobian matrix J is an M×N matrix, where each component is expressed as:

[0054]

[0055] Step 3: Calculate the Jacobian matrix

[0056] The Jacobian matrix is ​​calculated using the adjoint reciprocity theorem:

[0057] Introducing accompanying electrical sources and accompanying magnetic source The following Maxwell's equations are obtained:

[0058]

[0059]

[0060] Considering the boundary conditions and substituting them into the above equation, we obtain the following equation:

[0061]

[0062] The above formula is the basic formula for solving the Jacobian matrix of the electromagnetic field. The Jacobian matrix is ​​calculated by selecting a suitable adjoint field source. If it is necessary to obtain the E at the observation point r0... x Jacobi style, making and Then we get:

[0063]

[0064] The above formula is E x The formula for calculating the Jacobian of a component is used to obtain... The electric field value E can be obtained by first solving for the original electromagnetic field, and then by calculating the electric field value E under the condition of the accompanying source. * Finally, the dot product E * E is sufficient.

[0065] Step 4: Joint Inversion Imaging

[0066] The fitted function in the objective function is expanded to:

[0067] ||W(dF(m))|| 2 =||W1(d1-F1(m))|| 2 +||W2(d2-F2(m))|| 2

[0068]

[0069] Assuming the desired fitting error is achieved by increasing the data weighting coefficients α1 and α2, where the values ​​of α1 and α2 are related to the amount of data in each subset, the normalized joint fitting function is as follows:

[0070]

[0071] In the formula It is a data fitting balance weight that normalizes the fitting function for each data subset so that the smaller data subset has the same impact on the overall error as the larger subset.

[0072] The frequency-domain ocean controllable source electromagnetic and magnetotelluric data joint imaging method of the present invention has the following beneficial effects:

[0073] 1. Compared with traditional electromagnetic data imaging methods, the technology of this invention can make up for the shortcomings of single dataset inversion imaging, and the imaging results have higher accuracy and higher resolution.

[0074] 2. Compared with traditional electromagnetic data imaging methods, the present invention uses an unstructured dual-grid system to divide the model, with one set being a forward modeling grid and the other a reverse modeling grid. This improves the ability of marine controlled-source electromagnetic and magnetotelluric data imaging to process actual undulating terrain and complex media, and also offers good flexibility.

[0075] In summary, this invention overcomes the problems of insufficient imaging accuracy and low resolution in conventional marine controlled-source electromagnetic and magnetotelluric data imaging methods. This method can not only reconstruct actual undulating terrain and complex medium models, but also compensate for the shortcomings of inversion imaging using single datasets, significantly improving the accuracy and resolution of the imaging results. This has important practical value for guiding the inversion and interpretation of frequency-domain marine controlled-source electromagnetic and magnetotelluric data in complex regions. Attached Figure Description

[0076] Figure 1 This is an imaging flowchart of the frequency domain ocean controllable source electromagnetic and magnetotelluric data joint imaging method in an embodiment of the present invention.

[0077] Figure 2 This is a schematic diagram of mesh refinement.

[0078] Figure 3 This is the flowchart of the inversion algorithm.

[0079] Figure 4 This is an example of model inversion imaging. Detailed Implementation

[0080] The present invention will be further described in detail below with reference to the embodiments and accompanying drawings.

[0081] Example

[0082] The main steps in implementing this method include mesh refinement strategy, fast Occam's algorithm, Jacobian matrix calculation, and joint inversion imaging algorithm. The specific steps are as follows:

[0083] Step 1: Mesh Refinement Strategy. An unstructured dual-mesh system is adopted: one set is the inversion mesh, i.e., the coarse mesh, used for discretizing the parameters of the inversion model; the other set is the forward mesh, i.e., the fine mesh, used for discretizing the forward model to ensure the accuracy of the forward model, and is derived from the inversion mesh. Figure 2 a is a cross-sectional view of the initial mesh of a model. Figure 2 b is the cross-sectional view of its refined mesh.

[0084] The refinement strategy for the forward mesh in this invention is described below:

[0085] (1) Refine the cells near the receiving point, find the cell where the receiving point is located, and determine whether its volume is less than the set requirement for the volume of the receiving point cell. If it does not meet the requirement, halve its volume and refine it in the subsequent mesh refinement. Repeat the above steps until the volume of the cell meets the requirement.

[0086] (2) For electromagnetic data from a controlled marine source, in order to overcome the singularity of the emission source, the cells near the emission point need to be refined. Locate the cell where the emission point is located and determine whether its volume is smaller than the set requirement for the volume of the emission point cell. If it does not meet the requirement, its volume is halved and refined in the subsequent mesh refinement. Repeat the above steps until the volume of the cell meets the requirement.

[0087] (3) Iterate through each cell of the mesh to find the minimum distance from the cell to the receiver and the transmitter. Based on the distance and skin depth, calculate the maximum volume that the cell should meet. If the volume exceeds the maximum volume, its volume is halved and refined in the subsequent mesh refinement. Repeat the above steps until the volume of all cells meets the requirements.

[0088] (4) For inversion, in order to achieve unbiased transfer of inversion parameters to forward parameters, it is only necessary to mark each unit that needs to be inverted in the inversion grid as an independent region.

[0089] Step 2: Fast Occam's algorithm. The flowchart for the Fast Occam's inversion imaging algorithm is shown below. Figure 3 The algorithm constructs the inversion objective function using the Tikhonov regularization method, and the functional U is expressed as:

[0090] U=μ||Rm|| 2 +||W(dF(m))|| 2

[0091] In the formula: m is an N-dimensional model parameter vector, typically the resistivity value; R is the roughness operator matrix; μ is the Lagrange multiplier, used to balance model roughness and data fitting error. When μ takes a large value, the inversion mainly searches for a smooth model; otherwise, it mainly searches for the minimum fitting error. W is the associated diagonal weighted matrix; d is the observed data vector; and F(m) is the forward response corresponding to model m.

[0092] Given an initial model m k The objective function is minimized using the following iterative method:

[0093]

[0094] The data correction vector is:

[0095]

[0096] The Jacobian matrix J is an M×N matrix, where each component is expressed as:

[0097]

[0098] Step 3: Jacobian matrix calculation. The Jacobian matrix is ​​calculated using the adjoint reciprocity theorem:

[0099] Introducing accompanying electrical sources and accompanying magnetic source The following Maxwell's equations are obtained:

[0100]

[0101]

[0102] Considering the boundary conditions and substituting them into the above equation, we obtain the following equation:

[0103]

[0104] The above formula is the basic formula for solving the Jacobian matrix of the electromagnetic field. The Jacobian matrix is ​​calculated by selecting a suitable adjoint field source. If it is necessary to obtain the E at the observation point r0... x Jacobi style, making and Then we get:

[0105]

[0106] The above formula is E x The formula for calculating the Jacobian of a component is used to obtain... The electric field value E can be obtained by first solving for the original electromagnetic field, and then by calculating the electric field value E under the condition of the accompanying source. * Finally, the dot product E *E is sufficient.

[0107] Step 4: Joint Inversion Imaging Algorithm. The fitting function in the objective function is extended to:

[0108] ||W(dF(m))|| 2 =||W1(d1-F1(m))|| 2 +||W2(d2-F2(m))|| 2

[0109]

[0110] Assuming the desired fitting error is achieved by increasing the data weighting coefficients α1 and α2, where the values ​​of α1 and α2 are related to the amount of data in each subset, the normalized joint fitting function is as follows:

[0111]

[0112] In the formula It is a data fitting balance weight that normalizes the fitting function for each data subset so that the smaller data subset has the same impact on the overall error as the larger subset.

[0113] To illustrate the imaging effect of the method of the present invention, this embodiment takes a complex structural model as an example and performs controlled-source electromagnetic data inversion, magnetotelluric data inversion, and joint inversion imaging of controlled-source electromagnetic and magnetotelluric data respectively. Figure 4 a represents the inversion imaging result of a single magnetotelluric data set. Figure 4 b represents the inversion imaging result of electromagnetic data from a single controllable source. Figure 4 c represents the joint inversion imaging result of controlled-source electromagnetic and magnetotelluric data. From Figure 4 As can be seen from the data, the magnetotelluric data inversion imaging results show a more obvious reflection of low-resistivity bodies, but a poorer recovery of high-resistivity bodies. In contrast, the controlled-source electromagnetic data inversion imaging results recover the shallow high-resistivity anomalies and the undulations of the bedrock better, but the geometric shape and resistivity of low-resistivity anomalies are poorly recovered. The joint inversion imaging results accurately reflect the resistivity and location of the anomalies, without generating too many false anomalies, indicating that the joint inversion imaging results are reliable and have higher resolution. Figure 4 Figure d shows the variation curve of the root mean square (RMS) fit difference during the inversion imaging process. As can be seen from the figure, after multiple inversion iterations, the fit difference consistently converges to around 1, indicating that the inversion imaging results quickly and stably converge to the vicinity of the true model. This embodiment further illustrates that the method of the present invention has good imaging performance and can be applied to practical needs.

[0114] It will be apparent to those skilled in the art that the present invention is not limited to the details of the exemplary embodiments described above, and that the practical functions of the present invention can be realized in other specific forms without departing from the spirit or essential characteristics of the present invention.

Claims

1. A method for joint imaging of frequency-domain controlled-source electromagnetic and magnetotelluric data in the ocean, comprising the following steps: Step 1: Mesh Refinement An unstructured dual-grid system is used to discretize the forward and inverse model regions: one set is the inverse grid, i.e., the coarse grid, used to discretize the parameters of the inverse model; the other set is the forward grid, i.e., the fine grid, used to discretize the forward model to ensure the accuracy of the forward model, and is derived from the refinement of the inverse grid. Step 2: Fast Occam Inversion Imaging The inversion objective function is constructed using the Tikhonov regularization method, and the functional U is expressed as: U=μ||Rm|| 2 +||W(d-F(m))|| 2 In the formula: m is the N-dimensional model parameter vector, is the resistivity value, R is the roughness operator matrix, μ is the Lagrange multiplier, W is the associated diagonal weighted matrix, d is the observation data vector, and F(m) is the forward response corresponding to model m; Given an initial model m k The objective function is minimized using the following iterative method: The data correction vector is: The Jacobian matrix J is an M×N matrix, where each component is expressed as: Step 3: Calculate the Jacobian matrix The Jacobian matrix is ​​calculated using the adjoint reciprocity theorem: Introducing accompanying electrical sources and accompanying magnetic source The following Maxwell's equations are obtained: Considering the boundary conditions and substituting them into the above equation, we obtain the following equation: The above formula is the fundamental formula for solving the Jacobian matrix of the electromagnetic field. Let and Then we get: The above formula is E x The formula for calculating the Jacobian of the component involves first solving for the original electromagnetic field E, and then calculating the electric field value E in the case of the accompanying source. * Finally, the dot product E * ·E; Step 4: Joint Inversion Imaging The fitted function in the objective function is expanded to: ||W(d-F(m))|| 2 =||W1(d1-F1(m))|| 2 +||W2(d2-F2(m))|| 2 Assuming the desired fitting error is achieved by increasing the data weighting coefficients α1 and α2, where the values ​​of α1 and α2 are related to the amount of data in each subset, the normalized joint fitting function is as follows: In the formula It is a data fitting balance weight that normalizes the fitting function for each data subset so that the smaller data subset has the same impact on the overall error as the larger subset.

2. The frequency-domain ocean controlled-source electromagnetic and magnetotelluric data joint imaging method according to claim 1, characterized in that: The refinement steps of the forward mesh are as follows: Step 1.1 Refine the cells near the receiving point, find the cell where the receiving point is located, and determine whether its volume is smaller than the set requirement for the volume of the receiving point cell. If it does not meet the requirement, its volume is halved and then refined in the subsequent mesh refinement. Repeat the above steps until the volume of the cell meets the requirements; Step 1.2 If the data is from a controlled marine electromagnetic source, in order to overcome the singularity of the emission source point, the cells near the emission point need to be refined; find the cell where the emission point is located, and determine whether its volume is smaller than the set requirement for the volume of the emission point cell. If it does not meet the requirement, its volume is halved, and it is refined in the subsequent mesh refinement. Repeat the above steps until the volume of the cell meets the requirements; Step 1.3 Iterate through each cell of the mesh to find the minimum distance from the cell to the receiver and the transmitter. Based on the distance and skin depth, calculate the maximum volume that the cell should satisfy. If the volume exceeds the maximum volume, its volume is halved and refined in the subsequent mesh refinement. Repeat the above steps until the volume of all units meets the requirements; In step 1.4, during the inversion, each cell in the inversion grid that needs to be inverted is marked as an independent region.

Images