A method for combined magnetic-electric dual structure constrained airborne transient electromagnetic inversion

Abstract

The embodiment of the application discloses a kind of magnetic-electric double structure combined constraint's aviation transient electromagnetic inversion method, it is related to electromagnetic inversion imaging technical field, including: generating magnetic anomaly map based on aviation magnetic survey data, extracting magnetic boundary by edge detection algorithm, generating binary magnetism boundary image;Magnetic boundary coordinates in the magnetism boundary image are spatially registered with aviation transient electromagnetic inversion grid and projected to each depth horizon;Spatial structure framework matrix is constructed, and lateral structure constraint term is established;Based on ground transient electromagnetic sounding data inversion, shallow electrical section is obtained, and it is spatially registered with aviation transient electromagnetic inversion grid, shallow prior electrical constraint model is constructed and shallow prior electrical constraint term is determined;Magnetic-electric double constraint inversion objective function is constructed based on the above constraint term, iterative solution is carried out by inversion solving algorithm, and final inversion model is obtained after iteration termination.The application can improve the resolution of aviation electromagnetic method as a whole.

Description

Technical Field

[0001] This application relates to the field of electromagnetic inversion technology, and relates to, but is not limited to, an airborne transient electromagnetic inversion method with joint constraints of magnetic-electric dual structures. Background Technology

[0002] Airborne transient electromagnetic method (EEM) is a geophysical exploration method using aircraft as carriers. It has the advantages of wide coverage, no terrain limitations, and high efficiency, and is widely used in mineral resource exploration, hydrogeological surveys, and engineering surveys. However, the inversion interpretation of EEM faces two inherent challenges: insufficient lateral resolution and insufficient shallow resolution.

[0003] In existing technologies, airborne transient electromagnetic methods primarily detect anomalies by observing electrical differences in the subsurface medium. However, their lateral resolution is affected by factors such as power outage time, flight altitude, and transmit-receive distance, limiting their ability to delineate the horizontal boundaries of geological bodies and making it difficult to accurately determine the boundary locations and geometric morphologies of anomalies. Furthermore, since both the transmit and receive systems are located in the air, a large transmitting magnetic moment is typically required to acquire deep information. The turn-off time corresponding to a large current is slightly longer than that of a small current, causing the response signals of early time channels reflecting shallow information to overlap with the turn-off process, making separation difficult. Moreover, extracting the time channel information required for near-surface electrical structure using data processing techniques presents significant technical challenges, thus affecting the effectiveness of shallow exploration.

[0004] In contrast, airborne magnetics offers the advantage of high spatial resolution. Magnetic measurements are sensitive to the boundaries of magnetic bodies, enabling precise delineation of the spatial boundaries and geometry of underground magnetic bodies; however, its lateral resolution is relatively lower than that of airborne electromagnetics. While ground-based transient electromagnetic methods are less efficient than airborne methods, their receiving coils are close to the ground, providing high-resolution detection of shallow electrical structures and enabling the acquisition of clear shallow electrical profiles.

[0005] Therefore, a systematic structural constraint mechanism is urgently needed to effectively integrate the high lateral resolution of airborne magnetics with the high shallow resolution of ground transient electromagnetics, forming a joint constraint on airborne transient electromagnetic inversion and achieving a significant improvement in the overall detection capability of airborne electromagnetics. Summary of the Invention

[0006] This application provides an airborne transient electromagnetic inversion method with joint constraints of a magnetic-electric dual structure.

[0007] The technical solution of this application embodiment is implemented as follows: In a first aspect, embodiments of this application provide an airborne transient electromagnetic inversion method with joint constraints of magnetic and electric dual structures. The method includes: acquiring airborne magnetic survey data of an exploration area, preprocessing it to generate a magnetic anomaly map, extracting magnetic boundaries from the magnetic anomaly map using an improved Canny edge detection algorithm, and generating a binarized magnetic boundary image; spatially registering the magnetic boundary coordinates in the binarized magnetic boundary image with an airborne transient electromagnetic inversion grid and projecting them to various depth layers; constructing a spatial structure grid matrix and establishing lateral structural constraint terms from the magnetic boundaries to the electric spatial structure grid; acquiring ground transient electromagnetic bathymetry data of the exploration area and inverting to obtain shallow electric profiles; and then... The shallow electrical profile is spatially registered with the airborne transient electromagnetic inversion grid to construct a shallow prior electrical constraint model and determine the shallow prior electrical constraint terms. Based on the fitting difference term of the airborne transient electromagnetic data, the lateral structural constraint term, and the shallow prior electrical constraint term, a magnetic-electric dual-constraint inversion objective function is constructed. The magnetic-electric dual-constraint inversion objective function is iteratively solved using an inversion solution algorithm. In each iteration, the Jacobian matrix is ​​calculated and the model update amount is solved until the iteration termination condition is met, and the final inversion model is output. The final inversion model is compared with a control inversion model for lateral resolution evaluation, shallow resolution evaluation, and comprehensive evaluation to determine the dual-constraint inversion effect.

[0008] Optionally, the step of extracting the magnetic boundaries in the magnetic anomaly map using the improved Canny edge detection algorithm to generate a binary magnetic boundary image includes: performing convolution filtering on the magnetic anomaly map using a two-dimensional Gaussian kernel function to obtain a smoothed magnetic anomaly map; calculating the gradient magnitude and gradient direction of each pixel in the smoothed magnetic anomaly map; refining the edge width of the smoothed magnetic anomaly map based on the gradient magnitude and gradient direction using a non-maximum suppression algorithm to obtain a candidate edge image; and detecting and connecting edges using a dual thresholding method for the candidate edge image to generate a binary magnetic boundary image.

[0009] Optionally, in the spatial structure lattice matrix, grid nodes covered by magnetic boundaries are marked as boundary constraint active nodes, and the remaining grid nodes covered by non-magnetic boundaries are marked as boundary constraint inactive nodes; the establishment of the lateral structural constraint term from the magnetic boundary to the electric spatial structure lattice includes: calculating the boundary normal vector at each boundary constraint active node based on the spatial distribution of the boundary constraint active nodes in the airborne transient electromagnetic inversion grid; assigning boundary position weighting coefficients to each grid node in the horizontal two-dimensional space of the inversion region based on the spatial structure lattice matrix, wherein, at the boundary constraint active nodes, the boundary position weighting coefficients are set to a minimum value approaching 0; at the boundary constraint inactive nodes, the boundary position weighting coefficients are set to a value approaching 1; and constructing the lateral structural constraint term based on the boundary normal vector and the boundary position weighting coefficient, the expression of the lateral structural constraint term being represented by the following formula: ; In the formula, The model parameters to be inverted are represented as resistivity or conductivity. Ω represents the horizontal structural constraint term; Ω represents the horizontal two-dimensional space of the inversion region. Indicates the weighting coefficients for boundary positions; and Indicates model parameters in direction and First-order partial derivative in the direction; Represents the boundary normal vector. Represents the boundary quantity in the x-direction. This represents the boundary quantity in the y-direction.

[0010] Optionally, the shallow electrical profile has a depth of less than or equal to 150 m, and includes information on the electrical stratification structure and resistivity spatial distribution of the shallow strata; the expression for the shallow prior electrical constraint term is represented by the following formula: ; In the formula, This represents shallow prior electrical constraint terms; Indicates the parameters of the model to be inverted; This represents a shallow prior electrical constraint model; Indicates the inversion of three-dimensional space; (z) represents the prior constraint weighting coefficient, and the expression for the prior constraint weighting coefficient function is given by the following formula: ; In the formula, Represents the weighting coefficients of prior constraints; Represents depth coordinates; Indicates the initial weighting coefficients at the Earth's surface; This represents the attenuation constant.

[0011] Optionally, the objective function of the magneto-electric dual-constraint inversion is expressed by the following equation: ; In the formula, Represent the objective function for the magneto-electric dual-constraint inversion; This represents the fitting difference term in airborne transient electromagnetic data. , For airborne transient electromagnetic observation data, Indicates the parameters of the model to be inverted. This refers to the forward modeling response data obtained by performing forward modeling calculations based on the current parameters of the model to be inverted. Represents the data weight matrix; This represents a lateral structure constraint term, used to force the inversion results to remain spatially smooth in non-magnetic boundary regions, and to allow jumps in electrical parameters at magnetic boundary locations. , This represents a spatial structure constraint operator constructed based on aerospace magnetic boundaries; This represents a shallow prior electrical constraint term, used to force the shallow electrical properties of the inversion model to converge towards the ground-based transient electromagnetic prior model. , This represents a depth-weighted matrix, where the weights are constrained to decrease as depth increases. This represents a shallow prior electrical constraint model; This represents the first regularization parameter, used to control the relative weights of spatial structure constraints; This represents the second regularization parameter, used to control the relative weights of shallow prior constraints.

[0012] Optionally, the step of calculating the Jacobian matrix and solving for the model update amount until the iteration termination condition is met, and outputting the final inversion model, includes: approximating the Jacobian matrix using the finite difference method, wherein the calculation formula for each element in the Jacobian matrix is ​​expressed by the following formula: ; In the formula, Represents the Jacobian matrix in the form of the first... Line number Column elements; Indicates the first Forward modeling response function of airborne transient electromagnetic observation data; Indicates the first One model parameter; The perturbation step size represents the model parameters; singular value decomposition is performed on the Jacobian matrix, and the decomposition formula is expressed by the following equation: ; In the formula, Represents the Jacobian matrix; Represents the feature vector matrix of the data; Represents the model parameter eigenvector matrix; Indicates the transpose operation; Let represent a singular value diagonal matrix, where the singular values ​​satisfy ... ; Represents the maximum singular value. Indicates the second singular value. Represents the third singular value. Let represent the smallest non-zero singular value; set the basic damping factor, and calculate the model update amount for the current iteration step. The formula for calculating the model update amount is expressed as follows: ; In the formula, Indicates the amount of model updates; Represents the model parameter eigenvector matrix; Represents the identity matrix; Represents the data residual vector; Let represent the basic damping factor, which is adaptively adjusted during the iteration process, gradually decreasing from 0.1 to 0.01. The formula for calculating the basic damping factor is expressed as follows: ; In the formula, Indicates the first One basic damping factor; Indicates the basic damping factor; Represents singular value components The normalized coefficient, The model parameters are updated based on the model update amount, and it is determined whether any iteration termination condition is met. If the condition is met, the iteration is stopped, and the final inversion model is output.

[0013] Optionally, the iteration termination conditions include a root mean square fitting error of less than 1% and a model change between two adjacent iterations of less than a preset threshold, wherein the preset threshold is less than 1 × 10⁻⁶. -4 The current iteration count has reached the maximum iteration count, which is 200.

[0014] Secondly, embodiments of this application provide an electronic device, including a memory and a processor. The memory stores a computer program that can run on the processor. When the processor executes the program, it implements the steps in the above-described method for airborne transient electromagnetic inversion with combined magnetic-electric dual-structure constraints.

[0015] Thirdly, embodiments of this application provide a computer-readable storage medium storing a computer program thereon, which, when executed by a processor, implements the steps in the above-described airborne transient electromagnetic inversion method with combined magnetic-electric dual-structure constraints.

[0016] The beneficial effects of the technical solutions provided in this application include at least the following: This application provides an airborne transient electromagnetic inversion method with joint constraints of magnetic and electric dual structures. It acquires airborne magnetic survey data from the exploration area, preprocesses it to generate a magnetic anomaly map, and extracts magnetic boundaries from the magnetic anomaly map using an improved Canny edge detection algorithm to generate a binarized magnetic boundary image. The coordinates of the magnetic boundaries in the binarized magnetic boundary image are spatially registered with the airborne transient electromagnetic inversion grid and projected to each depth layer. A spatial structure grid matrix is ​​constructed, and a lateral structural constraint term is established from the magnetic boundaries to the electric spatial structure grid. Leveraging the high spatial resolution of airborne magnetic methods, the structural grid is constructed by extracting the magnetic boundaries and transformed into a lateral structural constraint for airborne transient electromagnetic inversion. This constraint effectively solves the inherent defect of insufficient vertical resolution in airborne electromagnetic methods and significantly improves the accuracy of longitudinal boundary positioning. Surface transient electromagnetic (TEM) bathymetry data of the exploration area were acquired, and shallow electrical profiles were obtained through inversion. The shallow electrical profiles were spatially registered with the airborne TEM inversion grid to construct a shallow prior electrical constraint model and determine the shallow prior electrical constraint terms. By employing an exponential decay weighting strategy, the surface electrical constraints were focused on the shallow layer, avoiding excessive extrapolation of shallow information to the deeper layers. This constraint effectively compensated for the insufficient resolution of shallow signals caused by airborne TEM at flight altitude. Based on a weighted combination of the fitting difference term, lateral structure constraint term, and shallow prior electrical constraint term from the airborne TEM data, a magnetic-electric dual-constraint inversion objective function was constructed. By setting regularization parameters, the spatial structure constraint and the shallow prior constraint were made independent. The system can be flexibly configured based on actual data availability. When only airborne magnetic data is available, spatial structure constraints can be used alone; when only ground electrical resistivity data is available, shallow constraints can be used alone, demonstrating good adaptability and scalability. An inversion algorithm iteratively solves the objective function of the magnetic-electric dual-constraint inversion. In each iteration, the Jacobian matrix is ​​calculated and the model update is determined until the iteration termination condition is met, outputting the final inversion model. During the iteration process, a differentiated damping strategy adaptively adjusts the basic damping factor, prioritizing the recovery of the model's macroscopic structure in the early stages of the inversion and gradually introducing fine model details in the later stages. This ensures inversion stability, significantly reduces inversion ambiguity, and improves the accuracy of lateral boundary positioning and shallow electrical property recovery. The final inversion model is compared with a control inversion model for lateral resolution evaluation, shallow resolution evaluation, and comprehensive evaluation to determine the effectiveness of the dual-constraint inversion.This application systematically integrates the high lateral resolution of airborne magnetics with the high shallow resolution of ground transient electromagnetics to form an inversion framework of "magnetic-electric dual-structure joint constraint". Unlike existing technologies that use the same type of data, this application utilizes the inherent complementary advantages of different geophysical methods in terms of resolution. Specifically, it extracts high-precision magnetic boundaries from airborne magnetic data to construct a spatial structural lattice, constraining the spatial geometric boundaries of airborne transient electromagnetic inversion and compensating for the insufficient spatial resolution of airborne electromagnetics. Furthermore, it integrates high-resolution shallow electrical profiles of ground transient electromagnetics to construct a shallow prior model, compensating for the insufficient shallow resolution of airborne electromagnetics, and achieving synergistic enhancement of airborne transient electromagnetic inversion in both lateral and vertical dimensions as well as in the shallow dimension. Attached Figure Description

[0017] To more clearly illustrate the technical solutions in the embodiments of this application, the accompanying drawings used in the description of the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of this application. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort, wherein: Figure 1 A flowchart illustrating an airborne transient electromagnetic inversion method with joint constraints of magnetic and electric dual structures, provided for an embodiment of this application; Figure 2 A schematic diagram of high lateral resolution boundary extraction and structural lattice construction using airborne magnetic methods provided in this application embodiment; Figure 3 This is a schematic diagram of the hardware entity of an electronic device provided in an embodiment of this application. Detailed Implementation

[0018] To make the objectives, technical solutions, and advantages of the embodiments of this application clearer, the technical solutions of the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of this application, not all embodiments. The following embodiments are used to illustrate this application, but are not intended to limit the scope of this application. Based on the embodiments in this application, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this application.

[0019] In the following description, references are made to “some embodiments,” which describe a subset of all possible embodiments. However, it is understood that “some embodiments” may be the same subset or different subsets of all possible embodiments and may be combined with each other without conflict.

[0020] It should be noted that the terms "first, second, and third" used in the embodiments of this application are merely to distinguish similar objects and do not represent a specific ordering of objects. It is understood that "first, second, and third" can be interchanged in a specific order or sequence where permitted, so that the embodiments of this application described herein can be implemented in an order other than that illustrated or described herein.

[0021] It will be understood by those skilled in the art that, unless otherwise defined, all terms used herein (including technical and scientific terms) have the same meaning as commonly understood by one of ordinary skill in the art to which the embodiments of this application pertain. It should also be understood that terms such as those defined in general dictionaries should be understood to have a meaning consistent with their meaning in the context of the prior art, and should not be interpreted in an idealized or overly formal sense unless specifically defined as herein.

[0022] The embodiments of this application will be further described below with reference to the accompanying drawings.

[0023] In view of the current problems in the study of airborne transient electromagnetic inversion in the field of electromagnetic inversion imaging technology, this application provides an airborne transient electromagnetic inversion method with joint constraints of magnetic-electric dual structures.

[0024] The technical solution of this application is described below, starting with the method embodiments.

[0025] Please refer to Figure 1 It shows a flowchart of an airborne transient electromagnetic inversion method with joint magnetic-electric dual-structure constraints provided in an embodiment of this application, as follows: Figure 1 As shown, the method includes at least the following steps S110 to S160.

[0026] Step S110: Obtain airborne magnetic survey data of the exploration area, generate a magnetic anomaly map after preprocessing, and extract the magnetic boundaries in the magnetic anomaly map using an improved Canny edge detection algorithm to generate a binarized magnetic boundary image.

[0027] In this embodiment, airborne magnetic survey data of the exploration area is acquired. This data includes various interfering components such as diurnal variation, flight altitude variation, and normal geomagnetic field background. To extract effective magnetic anomaly information reflecting the distribution of underground magnetic bodies, the airborne magnetic survey data undergoes diurnal variation correction, altitude correction, normal field correction, and magnetic anomaly extraction sequentially. Specifically, diurnal variation correction utilizes the diurnal variation curve of the geomagnetic field recorded by diurnal observation stations to remove the diurnal influence from the airborne magnetic survey data, eliminating time-varying interference from the external magnetic field caused by solar activity. Altitude correction unifies airborne magnetic survey data from different flight altitudes to the same reference surface to eliminate magnetic field variations caused by flight altitude changes. Normal field correction uses an IGRF model or geomagnetic field model to remove the Earth's main magnetic field background from the airborne magnetic survey data, highlighting magnetic anomalies caused by local geological anomalies. Magnetic anomaly extraction uses filtering or polynomial fitting methods to separate the regional field from the local field. After the above preprocessing is completed, a magnetic anomaly map reflecting the spatial distribution of underground magnetic bodies in the exploration area is obtained.

[0028] In this embodiment, the magnetic boundaries in the magnetic anomaly map are extracted using an improved Canny edge detection algorithm to generate a binary magnetic boundary image. Specifically, the magnetic anomaly map is smoothed by convolution filtering using a two-dimensional Gaussian kernel function to obtain a smoothed magnetic anomaly map. The expression for this two-dimensional Gaussian kernel function is as follows: ; In the formula, Represents a two-dimensional Gaussian kernel function; Represents the pixel coordinates of the image. These are the horizontal coordinates in the image coordinate system. These are the vertical coordinates in the image coordinate system; The Gaussian filter parameter is used to control the smoothing scale by adjusting its value, thereby suppressing the inherent high-frequency measurement noise in the magnetic anomaly map while preserving the main geometric features of the magnetic boundaries. In this embodiment, the Gaussian filter parameter is set to 1.5. In practical applications, it can be flexibly adjusted according to the noise level of the data in the magnetic anomaly map. When the noise is strong, the Gaussian filter parameter value is increased to enhance the smoothing effect; when the boundary detail requirement is high, the Gaussian filter parameter value is decreased to preserve the fine structure. Through the above Gaussian filter convolution operation, each pixel in the magnetic anomaly map is subjected to weighted smoothing processing to obtain a smoothed magnetic anomaly map.

[0029] Further, the gradient magnitude and gradient direction of each pixel in the smoothed magnetic anomaly map are calculated. Specifically, for each pixel in the smoothed magnetic anomaly map, the first-order partial derivatives of that pixel in the X and Y directions are calculated to obtain the gradient magnitude and gradient direction. The formulas for calculating the gradient magnitude and gradient direction are: ; ; In the formula, M represents the smooth magnetic anomaly map; This represents the gradient magnitude, used to characterize the severity of changes in magnetic anomalies. Indicates the gradient direction, used to characterize the direction of the fastest change in magnetic anomalies; This represents the rate of change of the smooth magnetic anomaly map in the x-direction; This represents the rate of change of the smooth magnetic anomaly map in the y-direction.

[0030] Furthermore, based on the gradient magnitude and gradient direction, the edge width of the smooth magnetic anomaly map is precisely refined using a non-maximum suppression algorithm to obtain candidate edge images. Specifically, all pixels in the smooth magnetic anomaly map are traversed along the gradient direction. For each pixel, it is determined whether the gradient magnitude of that point is a local maximum along the gradient direction; if so, the gradient magnitude of that point is retained; otherwise, the gradient magnitude of that point is set to zero.

[0031] Furthermore, for the candidate edge image, a dual-threshold method is used to detect and connect edges, generating a binary magnetic boundary image. Specifically, a high threshold and a low threshold are set, and the gradient magnitude of each pixel in the candidate edge image is compared with the high threshold and the low threshold. Pixels with gradient magnitudes higher than the high threshold are identified as strong edges and retained; pixels with gradient magnitudes lower than the low threshold are discarded; for weak edge pixels with gradient magnitudes between the low and high thresholds, they are retained only if they are spatially connected to strong edges, otherwise they are discarded. Through the above dual-threshold hysteresis connection strategy, a continuous and complete magnetic boundary is formed, ultimately generating a binary magnetic boundary image. In this embodiment, the high threshold is set to the magnitude corresponding to the 70% to 90% cumulative percentile of the gradient magnitude histogram, and the low threshold is set to 30% to 50% of the high threshold. The specific values ​​can be adaptively determined according to the gradient magnitude distribution characteristics of the magnetic anomaly map. Due to the high spatial resolution of airborne magnetics, the extracted magnetic boundary can reflect the spatial distribution range of underground geological bodies, providing high-quality spatial structural constraints for airborne transient electromagnetic inversion.

[0032] Step S120: Spatial registration of the magnetic boundary coordinates in the binarized magnetic boundary image with the airborne transient electromagnetic inversion grid, and projection onto each depth layer; constructing a spatial structure lattice matrix, and establishing lateral structural constraints from the magnetic boundary to the electric spatial structure lattice.

[0033] In this embodiment, the magnetic boundary coordinates in the binarized magnetic boundary image are spatially registered with the airborne transient electromagnetic inversion grid. Specifically, the spatial coordinates of each magnetic boundary point in the binarized magnetic boundary image are obtained, and the spatial range of the airborne transient electromagnetic inversion grid is obtained. This three-dimensional inversion grid spatial range is defined as the X-direction. , Y direction Z direction An affine transformation method is used to accurately map the magnetic boundary coordinates from the aerial survey image coordinate system to the Cartesian coordinate system of the airborne transient electromagnetic inversion grid, thereby achieving spatial registration between the magnetic boundary coordinates and the airborne transient electromagnetic inversion grid.

[0034] Furthermore, the registered magnetic boundary coordinates are projected onto each depth level of the airborne transient electromagnetic inversion grid. Specifically, for each depth level in the airborne transient electromagnetic inversion grid, the registered magnetic boundary coordinates are vertically projected onto the horizontal plane of that level, thus ensuring that the magnetic boundary maintains the same horizontal coordinates at each depth level.

[0035] In this embodiment, after projecting the magnetic boundary coordinates to each depth layer, a spatial structure lattice matrix is ​​constructed. Specifically, for each grid node in the horizontal two-dimensional space of the inversion region of the three-dimensional inversion grid, a value is assigned based on whether the node is covered by the magnetic boundary coordinates, thus generating the spatial structure lattice matrix. In the spatial structure lattice matrix, grid nodes covered by the magnetic boundary are marked as boundary constraint active nodes, and the remaining grid nodes not covered by the magnetic boundary are marked as boundary constraint inactive nodes. In this embodiment, the spatial extent of the airborne transient electromagnetic inversion grid is set to the X direction. , Y direction Z direction The inversion region is divided into a uniform grid with a grid cell size of [size missing]. ,common For each depth level, each grid cell is labeled. If a grid node is located on the magnetic boundary coordinates, it is marked as an active boundary constraint node and assigned a value of 1; otherwise, it is marked as an inactive boundary constraint node and assigned a value of 0. Each depth level independently constructs a spatial structure lattice matrix, and the magnetic boundary position maintains the same horizontal coordinates across all levels, thus achieving the vertical extension of the magnetic boundary into the deeper regions. Please refer to [reference needed]. Figure 2 It shows a schematic diagram of high lateral resolution boundary extraction and structural lattice construction by airborne magnetic method provided in the embodiments of this application.

[0036] Furthermore, a lateral structural constraint term is established from the magnetic boundary to the electrical spatial structure lattice. The specific steps are as follows: First, based on the spatial distribution of the boundary constraint activation nodes in the airborne transient electromagnetic inversion grid, the boundary normal vector at each boundary constraint activation node is calculated. Specifically, based on the arrangement orientation of the boundary constraint activation nodes in the horizontal two-dimensional space of the inversion region, the tangent direction of the magnetic boundary at each boundary constraint activation node is determined, and then the normal vector perpendicular to this tangent direction is calculated. This normal vector is perpendicular to the boundary orientation and points to both sides of the boundary. For any grid node in the spatial structure lattice matrix marked as a boundary constraint activation node, the tangent direction of the boundary orientation at that point is determined by analyzing the spatial connection relationship between this node and its adjacent boundary constraint activation nodes, and then the unit normal vector perpendicular to this tangent direction is calculated.

[0037] Furthermore, based on the spatial structure lattice matrix, boundary position weighting coefficients are assigned to each grid node in the horizontal two-dimensional space of the inversion region. Specifically, the core concept of spatial structure constraints is that airborne magnetics is sensitive to spatial structure boundaries. The extracted magnetic boundaries can accurately constrain the spatial variation locations of electrical parameters in airborne transient electromagnetic inversion. At magnetic boundary locations, discontinuities in electrical parameters should be allowed, corresponding to the interfaces of different lithological units. In non-boundary regions, electrical parameters should maintain lateral smoothness. Based on the above concept, at boundary constraint activated nodes, the boundary position weighting coefficients are set to a minimum value approaching 0 to allow discontinuous jumps in the parameters of the model to be inverted at that location; at boundary constraint inactive nodes, the boundary position weighting coefficients are set to a value approaching 1 to force the parameters of the model to be inverted to maintain lateral spatial smoothness at that location. In this embodiment, the boundary position weighting coefficient at the boundary constraint activated nodes is 0.001, and the boundary position weighting coefficient at the boundary constraint inactive nodes is 1.0.

[0038] Finally, based on the boundary normal vector and the boundary position weighting coefficients, a lateral structural constraint term is constructed. The expression for this lateral structural constraint term is as follows: ; In the formula, The model parameters to be inverted are represented as resistivity or conductivity. Indicates lateral structural constraints; Represents the horizontal two-dimensional space of the inversion region; This represents the weighting coefficient at the boundary position. It takes a minimum value close to 0 at the active nodes of the boundary constraint to allow abrupt changes in model parameters, and a value close to 1 at the inactive nodes of the boundary constraint to force lateral smoothing. and This represents the first-order partial derivatives of the model parameters in the X and Y directions; Let n represent the boundary normal vector. xRepresents the boundary quantity in the x-direction. This represents the boundary quantity in the y-direction. This lateral structure constraint term is used to force the inversion results to remain spatially smooth in non-magnetic boundary regions and to allow jumps in electrical parameters at magnetic boundary locations, so as to use the high spatial resolution of airborne magnetism to compensate for the insufficient lateral resolution of airborne transient electromagnetics.

[0039] Step S130: Obtain ground transient electromagnetic sounding data of the exploration area and invert to obtain shallow electrical profile; spatially register the shallow electrical profile with the airborne transient electromagnetic inversion grid, construct a shallow a priori electrical constraint model and determine the shallow a priori electrical constraint terms.

[0040] In this embodiment, ground transient electromagnetic sounding data is acquired by deploying transmitting and receiving coils on the surface, with a spacing of 50m between the measuring points, covering the main geological units of the exploration area. The ground transient electromagnetic sounding data is processed using a high-resolution inversion method to obtain shallow electrical profiles. Specifically, the Occam inversion method or the lateral constraint inversion method is used to perform high-resolution inversion processing on the ground transient electromagnetic sounding data, focusing on shallow high-resolution electrical profiles with a depth of less than or equal to 150m. These shallow electrical profiles clearly reflect the electrical stratification structure of the shallow strata in the exploration area (including topsoil, weathered layer, bedrock top surface, etc.) and the spatial distribution of resistivity at each layer.

[0041] In this embodiment, the shallow electrical profile is spatially registered with the airborne transient electromagnetic inversion grid. Specifically, since the coordinates of the ground transient electromagnetic measurement points and the airborne transient electromagnetic inversion grid are not completely coincident, Kriging interpolation or inverse distance weighted interpolation is used to map the electrical data of the ground measurement point coordinate system to each node of the airborne transient electromagnetic inversion grid, ensuring complete spatial coordinate consistency. For example, Kriging interpolation is used to interpolate the electrical values ​​of each ground measurement point location in the shallow electrical profile to each grid node of the airborne transient electromagnetic inversion grid, completing the spatial registration.

[0042] Furthermore, a shallow prior electrical constraint model is constructed. Specifically, based on the registered electrical profile data, corresponding electrical reference values ​​are assigned to the shallow nodes of the airborne transient electromagnetic inversion grid, generating a shallow prior electrical constraint model. This shallow prior electrical constraint model assigns electrical reference values ​​obtained from the inversion of ground transient electromagnetic data to each grid node in the shallow layer, and background resistivity values ​​to each grid node in the deeper layer, such as... .

[0043] Furthermore, the shallow-layer prior electrical constraint terms are determined. Specifically, ground transient electromagnetic methods have high resolution in shallow layers, and ground electrical method constraints are mainly concentrated in the shallow layer. As the depth increases, the reliability of ground electrical method constraints decreases, and the weighting coefficients should be reduced accordingly to avoid over-extrapolating shallow information to deeper layers and causing misleading results. Based on this, the prior constraint weighting coefficients for the shallow-layer prior electrical constraint terms are first determined. These prior constraint weighting coefficients vary with depth. Their expression is: ; In the formula, Represents the weighting coefficients of prior constraints; Represents depth coordinates; This represents the initial weighting coefficients at the Earth's surface, indicating the confidence level of the transient electromagnetic data at the Earth's surface. This represents the decay constant, used to control how quickly the weighted coefficients of the prior constraints decay with depth. The range of values ​​is The effective detection depth of transient electromagnetic fields on the ground is determined. In this embodiment, the initial weighting coefficient at the ground surface is set to 1.0, and the attenuation constant is set to... The prior constraint weighting coefficients decrease exponentially with increasing depth, causing the depth... The prior constraint weighting coefficients at the surface decay to the initial weighting coefficients at the surface. The following is This allows the weight of shallow prior constraints for transient electromagnetic fields on the ground to decrease with increasing depth, thus avoiding the over-extension of shallow electrical information to deeper layers and causing inversion misleading results.

[0044] Finally, based on the shallow prior electrical constraint model and the prior constraint weighting coefficients, a shallow prior electrical constraint term is constructed. The expression for the shallow prior electrical constraint term is: ; In the formula, This represents shallow prior electrical constraint terms; Indicates the parameters of the model to be inverted; This represents a shallow prior electrical constraint model; Indicates the inversion of three-dimensional space; The weighting coefficients for the prior constraints decrease exponentially with increasing depth. This shallow prior electrical constraint term is used to force the shallow electrical properties of the inversion model to converge with the ground transient electromagnetic prior model, so as to use the shallow high-resolution detection capability of ground transient electromagnetic to compensate for the insufficient shallow resolution of airborne transient electromagnetic due to flight altitude.

[0045] Step S140: Based on the fitting difference term of the airborne transient electromagnetic data, the weighted combination of the lateral structural constraint term and the shallow prior electrical constraint term, a magnetic-electric dual-constraint inversion objective function is constructed.

[0046] In this embodiment, airborne transient electromagnetic observation data is acquired, and forward modeling response data is obtained based on the current parameters of the model to be inverted. Then, a fitting difference term for the airborne transient electromagnetic data is constructed. This fitting difference term measures the deviation between the airborne transient electromagnetic observation data and the forward modeling response data. The expression for the fitting difference term is as follows: ; In the formula, This represents the fitting difference term in airborne transient electromagnetic data; This represents the data weight matrix, which is set according to the noise level of the observation data in each time channel. The earlier time channels have lower noise levels and are assigned higher weights, while the later time channels have higher noise levels and are assigned lower weights. This represents transient electromagnetic observation data from airborne sources. Indicates the parameters of the model to be inverted; This represents the forward response data obtained by forward modeling based on the current parameters of the model to be inverted.

[0047] Simultaneously, a lateral structure constraint term is obtained. This term forces the inversion results to maintain spatial smoothness in non-magnetic boundary regions and allows jumps in electrical parameters at magnetic boundary locations. Its discretization matrix form is as follows: ; In the formula, Indicates lateral structural constraints; This represents a spatial structure constraint operator constructed based on aerospace magnetic boundaries. This spatial structure constraint operator is obtained by discretizing the transverse structure constraint term using the finite difference or finite volume method. This represents the parameters of the model to be inverted.

[0048] Simultaneously, a shallow prior electrical constraint term is obtained, which is used to force the shallow electrical properties of the inversion model to approximate the ground transient electromagnetic prior model. Its discretized matrix form is as follows: ; In the formula, This represents shallow prior electrical constraint terms; This represents the depth-weighted matrix, where the constraint weights decrease as the depth increases. It is obtained by depth-weighted discretization of the prior constraint weighting coefficients. Indicates the parameters of the model to be inverted; This represents a shallow prior electrical constraint model.

[0049] Finally, by setting a first regularization parameter and a second regularization parameter, a weighted combination of the fitting difference term, lateral structure constraint term, and shallow prior electrical constraint term of the airborne transient electromagnetic data is constructed to establish a magnetic-electric dual-constraint inversion objective function. The expression for this magnetic-electric dual-constraint inversion objective function is as follows: ; In the formula, Represent the objective function for the magneto-electric dual-constraint inversion; This represents the fitting difference term in airborne transient electromagnetic data; Indicates lateral structural constraints; This represents shallow prior electrical constraint terms; This represents the first regularization parameter, used to control the relative weights of the lateral structural constraint terms; The second regularization parameter controls the relative weight of shallow prior electrical constraints. The regularization parameter is selected using the L-curve criterion or generalized cross-validation to balance the contribution of data fitting and each constraint. In this embodiment, the optimal values ​​of the first and second regularization parameters are determined using the L-curve criterion. Specifically, by plotting logarithmic curves of the data residual norm versus the constraint norm under different combinations of regularization parameters, and selecting the parameter values ​​corresponding to the inflection points, L-curve analysis determines that the first regularization parameter is 0.8 and the second regularization parameter is 1.2.

[0050] The magnetic-electric dual-constraint inversion objective function determined in this application is used as the sole optimization objective for the inversion iterative solution, ensuring that the inversion process is simultaneously controlled by the synergistic driving force of the lateral structural constraint term and the shallow prior electrical constraint term. Specifically, the lateral structural constraint term utilizes the high spatial resolution of airborne magnetics to solve the problem of insufficient spatial resolution of airborne electromagnetics; while the shallow prior electrical constraint term utilizes the high shallow resolution of ground transient electromagnetics to solve the problem of insufficient shallow resolution of airborne electromagnetics.

[0051] Step S150: The magnetic-electric dual-constraint inversion objective function is iteratively solved using an inversion solution algorithm. In each iteration, the Jacobian matrix is ​​calculated and the model update amount is solved until the iteration termination condition is met, and the final inversion model is output.

[0052] In this embodiment, the singular value decomposition (SVD) algorithm is used to iteratively solve the objective function of the magnetic-electric dual-constraint inversion. However, the solution method is not limited to SVD, and other methods can also be applied. Specifically, firstly, the initial parameters for the inversion solution are set, and the initial model parameters are set. In this embodiment, a homogeneous half-space initial model is used, and the initial resistivity value is set to... The initial basic damping factor is set to 0.1, the minimum basic damping factor is set to 0.01, the maximum number of iterations is set to 200, and the convergence threshold is set to... Set the current iteration number. Then, the iterative solution process begins. In each iteration, airborne transient electromagnetic forward modeling is performed based on the model parameters of the current iteration step. After obtaining the forward response data, the data residual vector is calculated. Furthermore, the Jacobian matrix is ​​approximated using the finite difference method. The formula for calculating each element of this Jacobian matrix is ​​as follows: ; In the formula, Represents the Jacobian matrix in the form of the first... Line number Column elements; Indicates the first Forward modeling response function of airborne transient electromagnetic observation data; Indicates the first One model parameter; This represents the perturbation step size of the model parameters. In this embodiment, the Jacobian matrix is ​​calculated using the forward difference method, and the perturbation step size is set to the current model parameter values. .

[0053] Furthermore, singular value decomposition is performed on the Jacobian matrix. The formula for singular value decomposition is: ; In the formula, Represents the Jacobian matrix; Represents the feature vector matrix of the data; Represents the model parameter eigenvector matrix; Indicates the transpose operation; Let represent a singular value diagonal matrix, where the singular values ​​satisfy ... , Represents the maximum singular value. Indicates the second singular value. Represents the third singular value. It represents the smallest non-zero singular value.

[0054] Furthermore, after obtaining the singular value decomposition results, the basic damping factor is set and the model update amount for the current iteration step is calculated. The formula for calculating the model update amount is: ; In the formula, Indicates the amount of model updates; Represents the model parameter eigenvector matrix; Represents a singular value diagonal matrix; This represents the squares of the elements of a singular value diagonal matrix; Represents the identity matrix; Represents the data residual vector; The base damping factor is represented by the following formula: ; In the formula, Indicates the first One basic damping factor; Indicates the basic damping factor; Represents singular value components The normalized coefficient, The aforementioned differentiated damping strategy ensures that the model parameter changes corresponding to large singular values ​​are greater than those corresponding to small singular values. Specifically, the model parameter changes corresponding to large singular values ​​are larger (the main update direction during inversion), while the changes corresponding to small singular values ​​are effectively suppressed to guarantee inversion stability. During the iteration process, the basic damping factor is adaptively adjusted, gradually decreasing from an initial value of 0.1 to 0.01. This prioritizes restoring the macroscopic structure of the model in the early stages of inversion, while allowing more small singular values ​​to participate in the inversion in the later stages to gradually introduce finer model details.

[0055] After calculating the model update amount, update the model parameters. ,in These are the current model parameters. These are the updated model parameters; This represents the model update amount. It checks if any iteration termination condition is met. If so, it stops iterating and outputs the final inverted model; otherwise, it returns to perform the forward calculation again and continues iterating. The iteration termination condition includes three items. The first iteration termination condition is that the root mean square fitting error is less than... In this embodiment, the root mean square fitting error decreased to [value missing] in the 127th iteration. If the first termination condition is met, the iteration terminates and the final inversion model is output; the second termination condition is that the change in the model between two consecutive iterations is less than a preset threshold, which is less than... The third iteration termination condition is that the current iteration count reaches the maximum iteration count, which is 200.

[0056] Step S160: Perform lateral resolution evaluation, shallow resolution evaluation and comprehensive evaluation on the final inversion model and the control inversion model to determine the dual-constraint inversion effect.

[0057] In this embodiment, the comparative inversion model includes an unconstrained inversion model, a magnetically constrained inversion model only, and an electrically constrained inversion model only. The comparative inversion model is obtained by adjusting the first regularization parameter and the second regularization parameter. Specifically, they are set respectively. Obtain the unconstrained inversion model; set up the following respectively. Obtain the magnetically constrained inversion model only; set up the following respectively Obtain the electrical method-only constrained inversion model.

[0058] The final inversion model was compared and evaluated with three sets of control inversion models. Lateral resolution evaluation included boundary positioning error and boundary sharpness. Boundary positioning error was determined by calculating the average offset distance between the recovered horizontal boundary and the actual boundary; boundary sharpness was determined by calculating the degree of restoration of resistivity contrast on both sides of the boundary. Shallow resolution evaluation included shallow resistivity restoration accuracy and thin-layer identification capability. Shallow resistivity restoration accuracy was calculated by comparing the inversion model with the actual model or known geological information. The root mean square error was determined; the thin-layer identification capability was determined by evaluating the recovery effect on shallow thin layers or lenses with a thickness not exceeding 150 μm. The comprehensive evaluation included the overall model recovery accuracy and the reduction in inversion ambiguity. The overall model recovery accuracy was determined by calculating the correlation coefficient between the full-space inversion results and the actual model; the reduction in inversion ambiguity was determined by comparing the stability of the inversion results under different initial models. Finally, based on the various evaluation indicators, the dual-constraint synergistic enhancement effect of the final inversion model can be quantitatively determined, verifying the effectiveness of the airborne transient electromagnetic inversion method with joint magnetic-electric dual-structure constraints provided in this application.

[0059] In summary, the airborne transient electromagnetic inversion method with joint magnetic-electric dual-structure constraints provided in this application acquires airborne magnetic survey data of the exploration area, generates a magnetic anomaly map after preprocessing, and extracts the magnetic boundaries in the magnetic anomaly map using an improved Canny edge detection algorithm to generate a binarized magnetic boundary image. The coordinates of the magnetic boundaries in the binarized magnetic boundary image are spatially registered with the airborne transient electromagnetic inversion grid and projected to each depth layer. A spatial structure grid matrix is ​​constructed, and a lateral structural constraint term is established from the magnetic boundaries to the electrical spatial structure grid. Utilizing the high spatial resolution of airborne magnetic methods, the structural grid is constructed by extracting the magnetic boundaries, transforming it into a lateral structural constraint for airborne transient electromagnetic inversion. This constraint effectively solves the inherent defect of insufficient vertical resolution in airborne electromagnetic methods and significantly improves the accuracy of longitudinal boundary positioning. Surface transient electromagnetic (TEM) bathymetry data of the exploration area were acquired, and shallow electrical profiles were obtained through inversion. The shallow electrical profiles were spatially registered with the airborne TEM inversion grid to construct a shallow prior electrical constraint model and determine the shallow prior electrical constraint terms. By employing an exponential decay weighting strategy, the surface electrical constraints were focused on the shallow layer, avoiding excessive extrapolation of shallow information to the deeper layers. This constraint effectively compensated for the insufficient resolution of shallow signals caused by airborne TEM at flight altitude. Based on a weighted combination of the fitting difference term, lateral structure constraint term, and shallow prior electrical constraint term from the airborne TEM data, a magnetic-electric dual-constraint inversion objective function was constructed. By setting regularization parameters, the spatial structure constraint and the shallow prior constraint were made independent. The system can be flexibly configured based on actual data availability. When only airborne magnetic data is available, spatial structure constraints can be used alone; when only ground electrical resistivity data is available, shallow constraints can be used alone, demonstrating good adaptability and scalability. An inversion algorithm iteratively solves the objective function of the magnetic-electric dual-constraint inversion. In each iteration, the Jacobian matrix is ​​calculated and the model update is determined until the iteration termination condition is met, outputting the final inversion model. During the iteration process, a differentiated damping strategy adaptively adjusts the basic damping factor, prioritizing the recovery of the model's macroscopic structure in the early stages of the inversion and gradually introducing fine model details in the later stages. This ensures inversion stability, significantly reduces inversion ambiguity, and improves the accuracy of lateral boundary positioning and shallow electrical property recovery. The final inversion model is compared with a control inversion model for lateral resolution evaluation, shallow resolution evaluation, and comprehensive evaluation to determine the effectiveness of the dual-constraint inversion.This application systematically integrates the high lateral resolution of airborne magnetics with the high shallow resolution of ground transient electromagnetics to form an inversion framework of "magnetic-electric dual-structure joint constraint". Unlike existing technologies that use the same type of data, this application utilizes the inherent complementary advantages of different geophysical methods in terms of resolution. Specifically, it extracts high-precision magnetic boundaries from airborne magnetic data to construct a spatial structural lattice, constraining the spatial geometric boundaries of airborne transient electromagnetic inversion and compensating for the insufficient spatial resolution of airborne electromagnetics. Furthermore, it integrates the high-resolution shallow electrical profile of ground transient electromagnetics to construct a shallow prior model, compensating for the insufficient shallow resolution of airborne electromagnetics, and achieving synergistic enhancement of airborne transient electromagnetic inversion in both lateral and shallow dimensions.

[0060] It should be noted that, in the embodiments of this application, if the above-mentioned airborne transient electromagnetic inversion method with combined magnetic-electric dual-structure constraints is implemented as a software functional module and sold or used as an independent product, it can also be stored in a computer-readable storage medium. Based on this understanding, the technical solution of the embodiments of this application, or the part that contributes to the related technology, 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 an electronic device to execute all or part of the methods described in the various embodiments of this application. The aforementioned storage medium includes various media capable of storing program code, such as USB flash drives, portable hard drives, read-only memory (ROM), magnetic disks, or optical disks. Thus, the embodiments of this application are not limited to any specific hardware and software combination.

[0061] Correspondingly, embodiments of this application provide a computer-readable storage medium storing a computer program thereon. When executed by a processor, this computer program implements the steps in the airborne transient electromagnetic inversion method with joint constraints of a magnetic-electric dual structure as described in any of the above embodiments. Correspondingly, embodiments of this application also provide a computer program product. When executed by a processor of an electronic device, this computer program product is used to implement the steps in the airborne transient electromagnetic inversion method with joint constraints of a magnetic-electric dual structure as described in any of the above embodiments.

[0062] Based on the same technical concept, this application provides an electronic device for implementing the airborne transient electromagnetic inversion method with joint constraints of magnetic-electric dual structures described in the above method embodiments. Figure 3 This is a hardware entity diagram of an electronic device provided in an embodiment of this application, such as... Figure 3As shown, the electronic device 300 includes a memory 310 and a processor 320. The memory 310 stores a computer program that can run on the processor 320. When the processor 320 executes the program, it implements the steps in the airborne transient electromagnetic inversion method with joint constraints of magnetic-electric dual structure as described in any embodiment of this application.

[0063] The memory 310 is configured to store instructions and applications executable by the processor 320, and can also cache data to be processed or already processed by the processor 320 and various modules in the electronic device (e.g., image data, audio data, voice communication data and video communication data), which can be implemented by flash memory or random access memory (RAM).

[0064] The processor 320 executes the program to implement the steps of an airborne transient electromagnetic inversion method with joint constraints of a magnetic-electric dual structure, as described above. The processor 320 typically controls the overall operation of the electronic equipment 300.

[0065] The aforementioned processor can be at least one of the following: Application Specific Integrated Circuit (ASIC), Digital Signal Processor (DSP), Digital Signal Processing Device (DSPD), Programmable Logic Device (PLD), Field Programmable Gate Array (FPGA), Central Processing Unit (CPU), Controller, Microcontroller, and Microprocessor. It is understood that other electronic devices can also implement the functions of the aforementioned processor, and this application does not specifically limit the specific implementation.

[0066] The aforementioned computer storage media / memory can be read-only memory (ROM), programmable read-only memory (PROM), erasable programmable read-only memory (EPROM), electrically erasable programmable read-only memory (EEPROM), magnetic random access memory (FRAM), flash memory, magnetic surface memory, optical disc, or compact disc read-only memory (CD-ROM), etc.; or it can be various electronic devices that include one or any combination of the above-mentioned memories, such as mobile phones, computers, tablet devices, personal digital assistants, etc.

[0067] It should be noted that the descriptions of the storage medium and device embodiments above are similar to the descriptions of the method embodiments above, and have similar beneficial effects. For technical details not disclosed in the storage medium and device embodiments of this application, please refer to the descriptions of the method embodiments of this application for understanding.

[0068] It should be understood that the phrase "one embodiment" or "an embodiment" throughout the specification means that a specific feature, structure, or characteristic related to the embodiment is included in at least one embodiment of this application. Therefore, "in one embodiment" or "in an embodiment" appearing throughout the specification does not necessarily refer to the same embodiment. Furthermore, these specific features, structures, or characteristics can be combined in any suitable manner in one or more embodiments. It should be understood that in the various embodiments of this application, the sequence numbers of the above-described processes do not imply a sequential order of execution; the execution order of each process should be determined by its function and internal logic, and should not constitute any limitation on the implementation process of the embodiments of this application. The sequence numbers of the above-described embodiments are merely descriptive and do not represent the superiority or inferiority of the embodiments.

[0069] It should be noted that, in this document, the terms "comprising," "including," or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or apparatus that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such a process, method, article, or apparatus. Unless otherwise specified, an element defined by the phrase "comprising one..." does not exclude the presence of other identical elements in the process, method, article, or apparatus that includes that element.

[0070] In the several embodiments provided in this application, it should be understood that the disclosed devices and methods can be implemented in other ways. The device embodiments described above are merely illustrative. For example, the division of units is only a logical functional division, and in actual implementation, there may be other division methods, such as: multiple units or components can be combined, or integrated into another system, or some features can be ignored or not executed. In addition, the coupling, direct coupling, or communication connection between the various components shown or discussed can be through some interfaces, and the indirect coupling or communication connection between devices or units can be electrical, mechanical, or other forms.

[0071] The units described above as separate components may or may not be physically separate. The components shown as units may or may not be physical units. They may be located in one place or distributed across multiple network units. Some or all of the units may be selected to achieve the purpose of the embodiments of this application, depending on actual needs.

[0072] In addition, each functional unit in the various embodiments of this application can be integrated into one processing unit, or each unit can be a separate unit, or two or more units can be integrated into one unit; the integrated unit can be implemented in hardware or in the form of hardware plus software functional units.

[0073] Alternatively, if the integrated units described above are implemented as software functional modules and sold or used as independent products, they can also be stored in a computer-readable storage medium. Based on this understanding, the technical solutions of the embodiments of this application, or the parts that contribute to related technologies, 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 the device automatic test line to execute all or part of the methods described in the various embodiments of this application. The aforementioned storage medium includes various media capable of storing program code, such as mobile storage devices, ROMs, magnetic disks, or optical disks.

[0074] The methods disclosed in the several method embodiments provided in this application can be arbitrarily combined without conflict to obtain new method embodiments.

[0075] The features disclosed in the several method or device embodiments provided in this application can be arbitrarily combined without conflict to obtain new method or device embodiments.

[0076] The above description is merely an embodiment of this application, but the scope of protection of this application is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application. Therefore, the scope of protection of this application should be determined by the scope of the claims.

Claims

1. A method for airborne transient electromagnetic inversion with joint constraints of magnetic-electric dual structures, characterized in that, The method includes: Airborne magnetic survey data of the exploration area is acquired, and after preprocessing, a magnetic anomaly map is generated. The magnetic boundaries in the magnetic anomaly map are extracted by an improved Canny edge detection algorithm to generate a binary magnetic boundary image. The magnetic boundary coordinates in the binarized magnetic boundary image are spatially registered with the airborne transient electromagnetic inversion grid and projected to each depth layer; a spatial structure lattice matrix is ​​constructed, and a lateral structural constraint term is established from the magnetic boundary to the electric spatial structure lattice. Obtain surface transient electromagnetic bathymetry data of the exploration area and invert to obtain shallow electrical profiles; spatially register the shallow electrical profiles with the airborne transient electromagnetic inversion grid, construct a shallow prior electrical constraint model and determine the shallow prior electrical constraint terms; A magnetic-electric dual-constraint inversion objective function is constructed by weighting the fitting difference term of the airborne transient electromagnetic data, the lateral structural constraint term, and the shallow prior electrical constraint term. The magnetic-electric dual-constraint inversion objective function is iteratively solved by an inversion solution algorithm. In each iteration, the Jacobian matrix is ​​calculated and the model update amount is solved until the iteration termination condition is met, and the final inversion model is output. The final inversion model is compared with the control inversion model by performing lateral resolution evaluation, shallow resolution evaluation, and comprehensive evaluation to determine the dual-constraint inversion effect.

2. The method according to claim 1, characterized in that, The step of extracting the magnetic boundaries in the magnetic anomaly map using the improved Canny edge detection algorithm to generate a binarized magnetic boundary image includes: The magnetic anomaly map is smoothed by convolution filtering using a two-dimensional Gaussian kernel function to obtain a smoothed magnetic anomaly map. Calculate the gradient magnitude and gradient direction of each pixel in the smooth magnetic anomaly map; Based on the gradient magnitude and gradient direction, the edge width of the smooth magnetic anomaly map is refined using a non-maximum suppression algorithm to obtain a candidate edge image. For the candidate edge image, a dual thresholding method is used to detect and connect edges to generate a binarized magnetic boundary image.

3. The method according to claim 1, characterized in that, In the spatial structure lattice matrix, the grid nodes covered by the magnetic boundary are marked as boundary constraint active nodes, and the remaining grid nodes covered by the non-magnetic boundary are marked as boundary constraint inactive nodes. The lateral structural constraint terms for establishing the magnetic boundary to the electric spatial structure lattice include: Based on the spatial distribution of the boundary constraint activation nodes in the airborne transient electromagnetic inversion grid, calculate the boundary normal vector at each boundary constraint activation node; Based on the spatial structure lattice matrix, boundary position weighting coefficients are assigned to each grid node in the horizontal two-dimensional space of the inversion region. Specifically, at the boundary constraint activated nodes, the boundary position weighting coefficients are set to a minimum value close to 0; at the boundary constraint inactive nodes, the boundary position weighting coefficients are set to a value close to 1. Based on the boundary normal vector and the boundary position weighting coefficient, a lateral structural constraint term is constructed, and the expression of the lateral structural constraint term is represented by the following formula: ； In the formula, The model parameters to be inverted are represented as resistivity or conductivity. Ω represents the horizontal structural constraint term; Ω represents the horizontal two-dimensional space of the inversion region. Indicates the weighting coefficients for boundary positions; and Indicates model parameters in direction and First-order partial derivative in the direction; Represents the boundary normal vector. Represents the boundary quantity in the x-direction. This represents the boundary quantity in the y-direction.

4. The method according to claim 1, characterized in that, The shallow electrical profile has a depth of less than or equal to 150 m, and includes information on the electrical stratification structure and resistivity spatial distribution of the shallow strata; the expression for the shallow prior electrical constraint term is given by the following formula: ； In the formula, This represents shallow prior electrical constraint terms; Indicates the parameters of the model to be inverted; This represents a shallow prior electrical constraint model; Indicates the inversion of three-dimensional space; Let represent the weighting coefficients of the prior constraints. The expression for the prior constraint weighting coefficient function is given by the following formula: ； In the formula, Represents the weighting coefficients of prior constraints; Represents depth coordinates; Indicates the initial weighting coefficients at the Earth's surface; This represents the attenuation constant.

5. The method according to claim 1, characterized in that, The objective function for the magneto-electric dual-constraint inversion is expressed by the following equation: ； In the formula, Represent the objective function for the magneto-electric dual-constraint inversion; This represents the fitting difference term in airborne transient electromagnetic data. , For airborne transient electromagnetic observation data, Indicates the parameters of the model to be inverted. This refers to the forward modeling response data obtained by performing forward modeling calculations based on the current parameters of the model to be inverted. Represents the data weight matrix; This represents a lateral structure constraint term, used to force the inversion results to remain spatially smooth in non-magnetic boundary regions, and to allow jumps in electrical parameters at magnetic boundary locations. , This represents a spatial structure constraint operator constructed based on aerospace magnetic boundaries; This represents a shallow prior electrical constraint term, used to force the shallow electrical properties of the inversion model to converge towards the ground-based transient electromagnetic prior model. , This represents a depth-weighted matrix, where the weights are constrained to decrease as depth increases. This represents a shallow prior electrical constraint model; This represents the first regularization parameter, used to control the relative weights of spatial structure constraints; This represents the second regularization parameter, used to control the relative weights of shallow prior constraints.

6. The method according to claim 1, characterized in that, The process of calculating the Jacobian matrix and solving for the model update amount continues until the iteration termination condition is met, outputting the final inversion model, including: The Jacobian matrix is ​​approximated using the finite difference method, and the formula for calculating each element in the Jacobian matrix is ​​expressed by the following equation: ； In the formula, Represents the Jacobian matrix in the form of the first... Line 1 Column elements; Indicates the first Forward modeling response function of airborne transient electromagnetic observation data; Indicates the first One model parameter; This indicates the perturbation step size of the model parameters; The singular value decomposition of the Jacobian matrix is ​​expressed by the following formula: ； In the formula, Represents the Jacobian matrix; Represents the feature vector matrix of the data; Represents the model parameter eigenvector matrix; Indicates the transpose operation; Let represent a singular value diagonal matrix, where the singular values ​​satisfy ... ; Represents the maximum singular value. Indicates the second singular value. Represents the third singular value. Represents the smallest non-zero singular value; Define a basic damping factor and calculate the model update amount for the current iteration step. The formula for calculating the model update amount is expressed as follows: ； In the formula, Indicates the amount of model updates; Represents the identity matrix; Represents the data residual vector; Let represent the basic damping factor, which is adaptively adjusted during the iteration process, gradually decreasing from 0.1 to 0.

01. The formula for calculating the basic damping factor is expressed as follows: ； In the formula, Indicates the first One basic damping factor; Represents singular value components The normalized coefficient, ; The model parameters are updated based on the model update amount, and it is determined whether any iteration termination condition is met. If it is met, the iteration is stopped, and the final inversion model is output.

7. The method according to claim 6, characterized in that, The iteration termination conditions include a root mean square fitting error of less than 1% and a model change of less than a preset threshold between two consecutive iterations, wherein the preset threshold is less than 1 × 10⁻⁶. -4 The current iteration count has reached the maximum iteration count, which is 200.

8. An electronic device comprising a memory and a processor, the memory storing a computer program executable on the processor, characterized in that, When the processor executes the program, it implements the steps of the method according to any one of claims 1 to 7.

9. A computer-readable storage medium having a computer program stored thereon, characterized in that, When executed by a processor, the computer program implements the steps of the method according to any one of claims 1 to 7.

Images